JUL  -9 


Final  Report 
on 

Liquid-Propellant  Droplet  Combustion  and  Cluster  Behavior  at 
Supercritical  Conditions 

for 

AFOSR  Contract/Grant  F49620-98- 1  -0034 


Prepared  by 
Vigor  Yang 

Propulsion  Engineering  Research  Center 
and 

Department  of  Mechanical  Engineering 
The  Pennsylvania  State  University 
University  Park,  PA  16802 


Submitted  to: 

Air  Force  Office  of  Scientific  Research 
801  N.  Randolph  Street,  Room  732 
Arlington,  VA  22203-1977 


May  2001 


DISTRIBUTION  STATEMENT  A 

Approved  for  Public  Release 
Distribution  Unlimited 


20010731  039 


REPORT  DOCUMENTATION  PAGE 

_  AFRL-SR-BL-TR-Ol- 

Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions 

data  needed,  and  completing  and  reviewing  this  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  ^ 

this  burden  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-0188),  121  7 

4302  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comp  t  - 

valid  0MB  control  number.  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 

1.  REPORT  DATE  (DD-MM-YYYY) 
01/05/2001 

2.  REPORT  TYPE 

Final  Technical  Report 

3.  DATES  COVERED  (From  -  To) 

^  01/11/1997-10/31/2000 

4.  TITLE  AND  SUBTITLE 

5a.  CONTRACT  NUMBER 

Liquid-Propellant  Droplet  Combustion  and  Cluster  Behavior  at  Supercritical  Conditions 

5b.  GRANT  NUMBER 

F49620-98-1-0034 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

Vigor  Yang 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME{S)  AND  ADDRESS(ES) 

The  Pennsylvania  State  University 

104  Research  Building  East 

University  Park,  PA  16802 

8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 

9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Dr.  Mitat  A.  Birkan 

Air  Force  Office  of  Scientific  Research 

10.  SPONSOR/MONITOR’S  ACRONYM(S) 
AFOSR 

801  N.  Randolph  Street,  Room  732 
Arlington,  VA  22203-1977 

11.  SPONSOR/MONITOR’S  REPORT 
NUMBER(S) 

1 

12.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 


Approved  for  public  release;  distribution  unlimited. 


13.  SUPPLEMENTARY  NOTES 


14.  ABSTRACT 

A  systematic  investigation  of  supercritical  droplet  vaporization  and  cluster  behavior  has  been  conducted  based  on  the  complete 
conservation  equations  in  both  the  gas  and  liquid  phases.  The  research  work  addresses  a  variety  of  fundamental  issues  related  to  droplet 
vaporization  and  dynamics  at  realistic  conditions  typical  of  liquid>propellant  rocket  combustion  devices.  A  unified  treatment  of  real-fluid 
thermodynamics  has  been  developed  based  on  fundamental  theories.  Special  attention  was  given  to  the  thermodynamic  non-ideality  and 
transport  anomaly  in  the  transcritical  regime.  A  series  of  calculations  has  been  performed  to  examine  the  cluster  behavior  of  liquid  oxygen 
(LOX)  droplets  in  both  sub-  and  super-critical  hydrogen  environments.  Results  show  that  pressure  has  strong  effect  on  droplet  interactions, 
while  the  temperature  effect  is  relatively  minor  at  high  pressures.  The  hydrogen  density  plays  a  decisive  role  in  determining  droplet 
interactions  through  its  influence  on  the  temperature  and  mass  fraction  gradients  at  the  LOX  droplet  surface.  The  characteristics  of  LOX 
droplet  vaporization  in  forced-convective  environments  has  also  been  studied.  A  dimensionless  parameter  which  represents  the 

ratio  of  aerodynamic  and  viscous  forces,  is  found  to  be  the  major  factor  determining  the  droplet  deformation  under  supercritical  conditions. 
Results  of  droplet  lifetime  are  well  correlated  as  a  function  of  the  initial  droplet  Reynolds  number  and  pressure.  Finally,  the  interactions 
between  two  droplets  moving  in  tandem  in  supercritical  convective  environments  were  investigated  in  detail. 

15.  SUBJECT  TERMS 

supercritical  combustion,  liquid  rocket  engines,  droplet  vaporization  and  combustion 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION 

OF  ABSTRACT 

18.  NUMBER 
OF  PAGES 

19a.  NAME  OF  RESPONSIBLE  PERSON 
211 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

uu 

19b.  TELEPHONE  NUMBER  (include  area 
code) 

(814)  863-1502 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std.  Z39.18 


Abstract 


A  systematic  investigation  of  supercritical  droplet  vaporization  and  cluster 
behavior  has  been  conducted  based  on  the  complete  conservation  equations  in  both  the 
gas  and  liquid  phases.  The  research  work  addresses  a  variety  of  fundamental  issues 
related  to  droplet  vaporization  and  dynamics  at  realistic  conditions  typical  of  liquid- 
propellant  rocket  combustion  devices. 

A  unified  treatment  of  real-fluid  thermodynamics  has  been  developed  based  on 
fundamental  theories.  Special  attention  was  given  to  the  thermodynamic  non-ideality  and 
transport  anomaly  in  the  transcritical  regime.  A  modified  Soave-Redlich-Kwong  (SRK) 
equation  of  state  was  utilized  to  derive  all  the  thermodynamic  correlations,  which  were 
then  incorporated  into  the  numerical  scheme  to  enhance  numerical  efficiency  and 
robustness.  In  addition,  a  preconditioning  scheme  with  the  dual  time-stepping  integration 
technique  was  implemented  to  render  the  numerical  algorithm  capable  of  treating  fluid 
flows  at  low  Mach  numbers.  The  fluid  thermal  conductivity  and  viscosity  were  estimated 
by  means  of  the  extended  corresponding  states  one-fluid  principle,  along  with  the 
Benedict-Webb-Rubin  (BWR)  equation  of  state.  The  pressure  effect  on  binary  mass 
diffusivity  is  corrected  by  the  Takahashi  method. 

A  series  of  calculations  has  been  performed  to  examine  the  cluster  behavior  of 
liquid  oxygen  (LOX)  droplets  in  both  sub-  and  super-critical  hydrogen  environments. 
Results  show  that  pressure  has  strong  effect  on  droplet  interactions,  while  the  temperature 
effect  is  relatively  minor  at  high  pressures.  The  hydrogen  density  plays  a  decisive  role  in 
determining  droplet  interactions  through  its  influence  on  the  temperature  and  mass 
fraction  gradients  at  the  LOX  droplet  surface.  In  the  later  stage  of  the  vaporization 


IV 


process,  droplet  vaporization  rate  is  dominated  by  the  oxygen  accumulation  at  the  droplet 
surface,  which  leads  to  low  thermal  conductivity  of  the  mixture  and  consequently  reduces 
the  heat  transfer  rate  into  the  liquid  droplet. 

Much  effort  was  expended  to  explore  LOX  droplet  vaporization  in  forced- 

convection  environments.  A  dimensionless  parameter  We^l^joh,  which  represents  the 
ratio  of  aerodynamic  and  viscous  forces,  is  found  to  be  the  major  factor  determining  the 
droplet  deformation  under  supercritical  conditions.  Results  of  droplet  lifetime  are  well 
correlated  as  a  function  of  the  initial  Reynolds  number  and  pressure.  A  linear  relationship 
is  generated  for  droplet  velocity,  and  the  droplet  trajectories  calculated  from  this  linear 
relationship  agree  well  with  the  results  directly  from  numerical  computation.  The  drag 
coefficient  in  the  supercritical  droplet  vaporization  process  was  found  to  decrease  with 
time  proportional  to  one-third  power  of  the  dimensionless  droplet  vaporization  rate, 
which  is  due  to  the  weak  Reynolds  number  reduction. 

The  interactions  with  two  droplets  moving  in  tandem  in  supercritical  convective 
environments  have  been  investigated  in  detail.  Results  indicate  that  droplet  dynamics 
exhibits  characteristics  distinct  from  that  of  an  isolated  droplet.  No  droplet  collision  was 
observed  for  all  the  cases  studied  herein.  A  forward  bag  break-up  of  the  leading  droplet, 
however,  was  found  when  the  two  droplets  are  initially  positioned  closely  with  a  H/R 
ratio  (the  ratio  of  the  initial  droplet  spacing  to  droplet  radius)  of  4  and  a  pressure  of  100 
atm.  The  presence  of  the  trailing  droplet  slightly  affects  the  lifetime  and  trajectory  of  the 
leading  droplets.  However,  the  lifetime  of  the  trailing  droplet  is  substantially  elongated. 
Forward  movement  was  observed  for  the  trailing  droplet  at  an  H/R  ratio  of  4  and  a 
pressure  of  lOOatm.  Increase  of  pressure  weakens  droplet  interactions. 
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Chapter  1 


INTRODUTION 


Spray  combustion  is  an  important  phenomenon  involved  in  various  liquid- 
propellant  combustion  devices,  including  diesel  engines,  gas-turbine  engines,  liquid- 
propellant  rocket  engines,  and  oil-fired  boilers.  During  spray  injection,  liquid  fuel  breaks 
up  into  small  droplets  through  atomization,  which  then  undergo  a  sequence  of 
vaporization,  mixing,  and  burning  processes.  The  droplet  vaporization  rate  has,  therefore, 
great  influence  on  the  engine  performance.  In  many  practical  liquid-fueled  combustors, 
the  chamber  pressures  and  temperatures  may  well  exceed  the  critical  conditions  of  the 
injected  fuels,  as  illustrated  in  Table  1.1.  Therefore,  droplet  vaporization  at  high  pressures 
has  received  considerable  attention  in  the  combustion  community. 


Table  1.1  Engine  Operation  Conditicms  and  Critical  Properties  of  Fuels 


Engines/Combustors 

Fuels 

Operation 

Conditions 

Critical  Properties  of 
Fuel 

Space  Shuttle 

Oxygen 

225  atm 

49.7  atm,  154.6K 

Main  Engine 

Hydrogen 

13  atm,  33.2  K 

Diesel  Engine 

Diesel  Oil 

>  50  atm 

~30  atm,  ~650K 

High-pressure  droplet  vaporization  exhibits  many  characteristics  distinct  from 
those  of  conventional  low-pressure  cases.  The  differences  lie  mainly  in  the  following 
areas.  First,  high-pressure  droplet  vaporization  is  transient  in  nature,  due  mainly  to  the 
rapid  droplet  regression  rate  and  the  transient  heating  of  the  droplet  interior.  When  the 
droplet  surface  reaches  the  critical  mixing  point,  the  sharp  distinction  between  the  liquid 
and  gaseous  phases  diminishes,  and  the  enthalpy  of  vaporization  decreases  down  to  zero, 
thereby  leading  to  a  high  droplet  vaporization  rate.  In  addition,  as  gaseous  properties  are 
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Strongly  modified  at  high  pressures,  the  diffusive  velocity  in  the  gaseous  region  becomes 
comparable  to  the  droplet  regressive  motion,  and  the  unsteadiness  is  further  enhanced. 
With  such  a  fast  vaporization  process,  a  uniform  temperature  distribution  inside  the 
droplet,  which  is  commonly  assumed  for  the  quasi-steady  case,  may  never  be  reached. 
The  transient  droplet  heating  may  persist  during  the  entire  droplet  lifetime,  especially 
when  supercritical  conditions  are  reached.  After  that,  there  exists  only  one  continuous 
fluid  medium  within  which  transient  heat  and  mass  transfer  take  place.  Second,  at  high 
pressures,  thermodynamic  non-ideality  becomes  an  important  factor,  which  means  the 
dissolution  of  ambient  gas  in  the  liquid  droplet  may  not  be  neglected,  and  the  partial 
molar  enthalpies  of  both  phases  must  be  used  to  calculate  the  enthalpy  of  vaporization. 
The  interfacial  thermodynamic  treatment  must  take  into  account  these  effects  by  means 
of  an  equation  of  state  for  real  fluids,  instead  of  the  ideal  equation  of  state.  Finally,  at 
high  pressures,  thermophysical  properties  become  functions  of  temperature,  pressure,  and 
compositions.  These  include  thermodynamic  properties,  such  as  enthalpy  and  the 
constant-pressure  heat  capacity,  and  transport  properties,  such  as  thermal  conductivity 
and  mass  diffusivity.  Furthermore,  the  property  anomaly  near  the  critical  state  becomes 
profound  and  must  be  treated  properly. 

1.1  Literature  Survey 

1.1.1  Droplet  Vaporization  in  Quiescent  Environments 

1.1. 1.1  Single  Droplet  Vaporization 

The  single  droplet  vaporization  phenomenon  has  been  studied  since  the  1950s. 
The  earliest  work  is  a  one-dimensional  low-pressure  analytical  solution  by  Godsave 
(1953),  which  was  based  on  a  few  assumptions,  including  spherical  symmetry,  quasi- 
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Steady  state,  constant  and  uniform  pressure,  etc.  Spalding  (1959)  studied  the  combustion 
of  a  single  fuel  droplet  at  supercritical  conditions,  with  an  assumption  that  the  droplet 
reaches  the  critical  point  immediately  after  being  injected  into  the  combustor.  The  droplet 
was  further  treated  as  a  point  source  of  mass.  The  analysis  was  later  improved  by  Rosner 
(1967)  to  include  a  finite  droplet  mass.  Both  studies  demonstrated  the  transient 
characteristics  of  droplet  burning  under  supercritical  conditions. 

Wieber  (1963)  investigated  a  single  droplet  vaporizing  in  a  convective 
environment  using  an  empirical  expression  of  the  Nusselt  number  to  calculate  the  heat 
transfer  rate  at  the  droplet  surface.  The  liquid  thermal  conductivity  was  taken  to  be 
infinity  by  assuming  rapid  internal  circulation  within  the  droplet,  which  also  implied  a 
uniform  distribution  of  temperature.  Both  heptane  and  oxygen  droplets  were  studied  with 
pressures  up  to  2000  psia.  Results  of  temperature  histories  concluded  that  the  droplet 
surface  temperature  eould  rise  continuously  and  finally  reach  the  critical  point. 

The  similar  problem  was  re-examined  by  a  number  of  researchers  (Brzustowski 
1965,  Chervinsky  1969,  Polymeropoulos  and  Peskin  1972,  Sanchez-Tarifa  et  al.,  1972)  in 
order  to  investigate  the  influences  of  convection,  density  variations,  and  finite-rate 
chemical  kinetics.  None  of  these  early  studies,  however,  provided  faithful  results, 
because  many  of  the  high-pressure  characteristics,  such  as  the  thermodynamic 
nonideality  and  transport  anomaly,  were  not  taken  into  account. 

A  systematic  investigation  of  high-pressure  phenomena  was  first  conducted  by 
Manrique  and  Borman  (1969).  In  their  work,  the  concepts  of  non-ideal  mixtures  and 
transport  property  anomaly  near  the  critical  mixing  point  were  introduced,  while  a  quasi¬ 
steady  vaporization  model  was  employed  to  simplify  the  analysis.  The  droplet  was 
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assumed  to  remain  spherical  symmetry  with  a  uniform  temperature  distribution  during 
the  entire  vaporization  process,  and  the  droplet  regression  rate  was  neglected.  At  the 
droplet  surface,  thermodynamic  equilibrium  was  considered,  with  the  dissolution  of  the 
ambient  gas  into  the  liquid  phase  confined  within  a  thin  layer  near  the  droplet  surface.  It 
was  demonstrated  that  the  difference  in  the  partial  molar  enthalpies  of  constitute  species 
between  the  gaseous  and  liquid  phases  must  be  included  in  the  energy  equation  to  replace 
the  latent  heat  of  vaporization  of  the  pure  liquid  species.  The  Redlich-Kwong  equation  of 
state  was  utilized  to  solve  the  two-component  thermodynamic  phase  equilibrium 
problem,  and  to  calculate  the  partial  molar  enthalpies.  The  carbon  dioxide-nitrogen 
system  was  studied  with  ambient  temperatures  and  pressures  ranging  from  500-1600K 
and  70-120  atm,  respectively.  The  mass  diffusivity  was  estimated  from  a  corresponding- 
state  chart,  and  the  thermal  conductivity  was  evaluated  by  using  the  Stiel  and  Thodos 
pure  component  correlation  (1964),  which  treats  the  mixture  as  a  hypothetical  pure 
component  with  pseudocritical  properties.  The  conclusion  was  that  the  thermodynamic 
non-ideality  effects,  such  as  the  dissolution  of  ambient  nitrogen  in  the  liquid  droplet  and 
application  of  the  partial  molar  enthalpies,  and  the  transport  property  anomaly  have 
significant  influences  at  high  pressures.  Large  errors  may  arise  without  considering  them 
properly.  For  a  given  ambient  temperature,  steady-state  vaporization  may  never  be 
attained  if  the  system  pressure  exceeds  a  threshold  value.  The  droplet  thus  vaporizes 
unsteadily  throughout  its  lifetime. 

Following  Manrique  and  Borman,  Matlosz  et  al.  (1972)  studied  n-hexane  droplet 
vaporization  in  a  stagnant  environment  of  nitrogen  and  argon  at  pressures  up  to  102  atm, 
where  the  thermodynamic  non-ideality  effects  on  phase  equilibrium  and  enthalpy  of 
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vaporization  were  included.  However,  all  the  other  thermophysical  properties  were  taken 
to  be  constant.  The  transient  behavior  of  droplet  vaporization,  such  as  surface  regression, 
was  taken  into  account,  but  the  temperature  inside  the  liquid  droplet  was  still  assumed  to 
be  uniform.  The  Redlich-Kwong  equation  of  state  was  chosen  to  treat  the  high-pressure 
non-ideal  mixture.  Numerical  results  of  the  time  histories  of  droplet  temperature  and 
radius  were  compared  with  experimental  measurements.  It  was  demonstrated  that  the 
non-ideal  behavior  of  the  gas  phase,  which  includes  the  dissolution  of  gaseous  species 
into  the  liquid  droplet,  is  important  at  high  pressures. 

To  clarify  the  systematic  errors  associated  with  the  quasi-steady  approximation  in 
the  droplet  vaporization  process  when  the  droplet  surface  approaches  its  critical  state, 
Rosner  and  Chang  (1973)  presented  exact  solutions  for  the  transient  vaporization  of  an 
isolated  spherical  droplet  in  a  quiescent  constant-property  environment.  Illustrative 
calculations  were  included  for  the  case  of  a  kerosene-like  (n-dodecane)  droplet 
evaporating  into  compression-heated  air,  similar  to  the  conditions  encountered  in  diesel 
engines.  The  non-ideality  effects  of  gas  dissolution  in  the  droplet  and  the  partial  molar 
enthalpy  were  not  considered,  but  it  was  indicated  that  the  effects  would  not  change  the 
conclusion  qualitatively.  With  pressures  up  to  2.96  times  of  the  critical  pressure,  the 
quasi-steady  theory  would  overestimate  the  droplet  lifetimes  by  47  percent,  and  the 
overestimation  would  reach  89  percent  with  the  presence  of  an  envelop  flame.  Results 
demonstrated  that  the  evaporation  process  for  an  isolated  fuel  droplet  near  its  critical 
condition  could  not  be  treated  as  quasi-steady. 

These  early  works  clearly  show  that  high-pressure  droplet  vaporization  is  an 
unsteady  process  during  which  a  droplet  is  continuously  heated  up,  and  the  critical  state 
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could  eventually  be  reached.  However,  when  a  thermodynamic  system  transits  from  its 
subcritical  to  supercritical  state,  the  singular  behavior  at  critical  state  remains  unclear. 
Umemura  (1986)  constructed  an  asymptotic  analysis  of  the  heat  and  mass  transfer  near 
the  gas-liquid  interface  at  the  critical  state  of  the  mixture.  He  concluded  that,  within  the 
range  of  validity  of  the  so-called  local  equilibrium  hypothesis  and  linear 
phenomenological  relations,  the  mutual  diffusion  coefficient  must  vanish  at  the  critical 
surface  to  ensure  a  finite  evaporation  rate  in  the  critical  condition. 

In  order  to  explore  the  natural  convection  effects  on  high-pressure  droplet 
vaporization,  Kadota  and  Hiroyasu  (1976)  investigated  a  single  n-heptane  droplet 
evaporating  in  high-pressure  and  high-temperature  nitrogen  environments.  In  this  study, 
thermodynamic  non-ideality  was  included,  and  thermophysical  properties  were  calculated 
by  fitted  polynomials  as  functions  of  temperature  and  pressure.  The  effect  of  natural 
convection  was  treated  using  empirical  relations  of  the  Nusselt  and  Sherwood  numbers. 
Jin  and  Borman  (1986)  applied  a  similar  model  to  study  the  vaporization  of 
multicomponent  fuel  droplets  of  pentane  and  octane  in  a  nitrogen  stream.  The  conclusion 
was  that  the  difference  of  the  evaporation  rates  between  the  two  species  decreases  with 
increasing  pressure,  and  the  potential  of  micro-explosion  is  inhibited  at  high  pressures. 

Hsieh  et  al.  (1991)  conducted  a  comprehensive  analysis  of  multicomponent  fuel 
droplet  vaporization  at  near  critical  conditions,  based  on  transient  conservation  equations 
with  a  thorough  treatment  of  property  variations.  Thermodynamic  non-ideality  effects  on 
liquid  dissolution  and  enthalpy  of  vaporization  were  modeled  through  the  use  of  a 
modified  Soave-Redlich-Kwong  (SRK)  equation  of  state,  which  was  proved  to  be  more 
accurate  than  the  conventional  Redlich-Kwong  equation  of  state.  A  moving  coordinate 
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system  was  employed  to  account  for  the  droplet  surface  regression.  Thermodynamic  and 
transport  properties,  such  as  constant-pressure  heat  capacity  and  thermal  conductivity, 
were  treated  as  functions  of  temperature  and  pressure  for  each  constituent  species,  and 
then  evaluated  for  a  mixture  using  appropriate  mixing  rules.  Pressure  effect  on  binary 
mass  diffusivity  was  estimated  using  the  Takahashi  (1974)  method.  The  analysis  was 
used  to  study  the  vaporization  process  of  single  droplets  containing  either  pure  n-pentane 
or  mixed  n-pentane/n-octane  in  gaseous  nitrogen  at  supercritical  conditions.  Calculations 
were  carried  out  through  either  the  end  of  droplet  lifetime  or  the  critical  point.  For 
multicomponent  fuel  systems,  the  droplet  vaporization  rate  increases  progressively  with 
ambient  pressure,  and  a  droplet  may  reach  its  critical  mixing  state  before  the  end  of  its 
lifetime  for  extreme  cases. 

Shuen  et  al.  (1992)  later  extended  the  work  by  Hsieh  et  al.  (1991)  to  study 
supercritical  droplet  combustion.  In  this  analysis,  complete  conservation  equations  along 
with  the  consideration  of  bulk  motion  in  the  liquid  phase  were  solved  using  the 
preconditioning  dual-timing  integration  technique  to  overcome  the  numerical  stiffness. 
Finite-rate  chemical  kinetics  was  accommodated  for  treating  the  gas-phase  combustion. 
Once  the  droplet  surface  reached  the  critical  mixing  state,  an  isotherm  line  at  the  critical 
mixing  temperature  was  employed  to  determine  the  gas-liquid  interface,  which  allowed 
the  model  to  compute  entire  droplet  vaporization  process.  A  n-paraffin  fuel  droplet 
burning  in  air  was  studied.  Results  for  large  droplets  (e.g.  1000  pm  in  diameter)  indicated 
that  at  subcritical  conditions,  the  vaporization  of  liquid  fuel  dominates  the  burning 
characteristics  of  the  droplet,  while  at  supercritical  conditions,  the  transient  diffusion  in 
the  gas  phase  controls  the  combustion  process.  The  combustion  lifetime  decreases  with 
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increasing  pressure,  and  reaches  a  minimum  near  the  critical  pressure  of  the  liquid  fuel. 
As  pressure  further  increases,  the  combustion  lifetime  increases  due  to  reduced  mass 
diffusivity  at  high  pressures,  which  is  in  qualitative  agreement  with  the  experimental 
observations  of  Faeth  et  al.  (1969)  and  Sato  et  al.  (1990). 

This  method  was  further  extended  by  Yang  et  al.  (1994)  to  treat  liquid  oxygen 
(LOX)  droplet  vaporization  in  a  hydrogen  environment.  A  modified  Soave-Redlich- 
Kwong  (SRK)  equation  of  state  was  utilized  with  special  attention  focused  on  the 
quantum-gas  behavior  of  hydrogen.  Significant  efforts  were  devoted  to  estimating  the 
thermophysical  properties  as  functions  of  pressure.  The  extended  corresponding  states 
model  of  Ely  and  Hanly  (1981,  1983)  was  applied  to  evaluate  thermal  conductivity,  heat 
capacity,  and  viscosity  of  real  fluid  mixtures.  Various  distinct  high-pressure  effects  on 
droplet  behavior  were  investigated  in  depth.  Results  indicated  that  droplet  lifetime  has 
strong  pressure  dependence,  and  can  be  well  correlated  with  the  square  of  the  initial 
droplet  diameter. 

Lafon  (1995)  studied  single  droplet  vaporization  in  a  gaseous  environment.  Both 
cases  of  liquid  oxygen  (LOX)  in  hydrogen  and  hydrocarbon  fuels  in  nitrogen  were 
investigated.  A  unified  treatment  of  thermodynamic  properties  and  numerical  algorithms 
for  real  fluid  mixtures  was  carried  out  based  on  fundamental  thermodynamic  theories. 
The  computational  efficiency  and  robustness  were  significantly  improved.  A  linear 
irreversible  theory  was  formulated  to  calculate  the  mass  vaporization  rate  at  the  droplet 
surface,  which  relates  the  vaporization  rate  to  the  difference  of  species  chemical  potential 
between  the  two  phases.  Extensive  parametric  studies  were  conducted,  and  results 
compared  well  with  experimental  measurements. 
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Haldenwang  et  al.  (1996)  also  investigated  the  vaporization  of  liquid  oxygen 
(LOX)  droplet  in  a  quiescent  hydrogen  environment.  The  transition  from  the  subcritical 
to  the  supercritical  vaporization  regime  was  examined.  After  the  droplet  reached  the 
supercritical  vaporization  regime,  the  droplet  surface  was  defined  based  on  the  critical 
mixing  mass  fraction  of  oxygen,  instead  of  the  critical  mixing  temperature  (Yang  et  al. 
1994).  A  minimum  of  droplet  lifetime  at  a  pressure  above  the  critical  pressure  of  the  fuel 
was  revealed. 

A  few  other  researchers  have  conducted  numerical  analyses  of  high-pressure 
droplet  vaporization  in  quiescent  environments,  notably  Litchford  and  Jeng  (1990),  Curtis 
and  Farrell  (1992),  Jia  and  Gogos  (1993),  Delplanque  and  Sirignano  (1993),  Umemura 
and  Shimada  (1996).  The  effect  of  ambient  pressure  oscillation  on  droplet  behavior  was 
also  investigated  (Lafon  et  al.  1995).  Two  recent  review  papers  are  available  (Yang  2000, 
and  Bellan  2000). 

The  research  works  mentioned  above  are  all  based  on  one  conventional 
assumption  that  thermodynamic  phase  equilibrium  prevails  at  the  droplet  surface,  which 
is  compatible  with  the  local  equilibrium  hypothesis  for  defining  local  thermodynamic 
variables,  such  as  temperature  and  density.  Recently,  an  attempt  was  made  to  investigate 
this  problem  using  non-equilibrium  theory.  Harstad  and  Bellan  (1998)  developed  a  model 
to  study  liquid  oxygen  (LOX)  droplet  vaporization  in  a  hydrogen  environment.  The 
conservation  equations  were  derived  based  on  the  fluctuation  theory  of  Keizer  (1987), 
which  account  for  both  Soret  and  Dufour  effects  in  the  calculation  of  heat  and  mass 
fluxes.  A  non-equilibrium  evaporation  law  was  incorporated  to  treat  the  mass 
vaporization  flux  at  the  droplet  surface,  which  calculates  the  flux  at  the  molecular  level. 
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In  this  work,  however,  the  non-equilibrium  formulation  played  no  role  in  the  entire 
computation,  because  only  supercritical  conditions  were  investigated.  The  same  model 
was  applied  to  study  n-heptane  droplet  vaporization  in  a  nitrogen  system  (Harstad  and 
Bellan,  1999).  Because  of  the  lack  of  methods  to  calculate  the  thermal  diffusion  factor  at 
high  pressures,  this  parameter  was  determined  with  arbitrariness.  After  the  droplet 
surface  reached  critical  condition,  it  was  defined  at  the  position  of  the  maximum  density 
gradient  in  order  to  be  compared  with  optical  measurements  (Sato  1993,  Nomura  1996, 
Morin  1999).  Reasonably  good  agreements  were  obtained  with  microgravity  experiments 
at  low  pressures  (from  0. 1-0.5  MPa),  but  obvious  differences  still  exist  at  high  pressures 
(1-2  MPa).  Because  these  microgravity  experiments  were  conducted  with  droplets 
suspended  on  a  fiber  in  parabolic  flights  or  drop  tower,  the  effects  of  buoyancy  and 
suspending  fiber  become  more  profound  at  high  pressures  (Vieille  et  al.  1996,  Morin 
1999).  It  is  not  clear  what  causes  the  differences  between  the  simulation  and  experimental 
results. 

Recently,  limited  success  has  also  been  achieved  for  exploring  supercritical 
droplet  vaporization  using  molecular  dynamics  (MD),  which  directly  simulates  the 
behavior  of  the  molecules  (Kaltz  et  al.  1998).  A  detailed  description  of  this  method  is 
beyond  the  scope  of  this  thesis.  Additional  information  about  the  current  progresses  in 
this  field  is  available  in  the  review  paper  by  Sirignano  and  Delplanque  (1999). 

1.1. 1.2  Droplet  Cluster  Behavior 

In  practical  spray  combustors,  there  is  no  evidence  that  a  droplet  behaves  in  an 
isolated  manner.  In  contrast,  droplets  normally  vaporize  and  burn  as  a  group  (McCreath 
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and  Chiger  1973).  Because  there  are  thousands  of  droplets  involved  in  the  problem, 
simplifications  are  required  in  numerical  analyses. 

Tishkoff  (1979)  studied  the  effect  of  droplet  interactions  by  modeling  the 
evaporation  of  a  droplet  inside  a  “bubble”.  The  “droplet-in-bubble”  concept  was 
developed  based  on  the  cellular  model  of  Zung  (1967),  and  considerable  improvement 
was  achieved.  The  droplets  were  assumed  to  be  monosized  and  distributed  uniformly  in  a 
group.  During  the  entire  vaporization  process,  all  the  droplets,  whether  on  the  boundary 
or  inside  the  cluster,  remained  a  spherical  shape  and  evaporated  in  a  quiescent 
environment.  The  difference  from  the  single  droplet  case  is  that  a  droplet  is  now 
evaporating  in  a  finite  space,  a  bubble.  The  heat  and  mass  transfer  between  the  droplet 
cluster  and  the  outside  environment  were  equally  distributed  to  each  bubble. 
Corresponding  to  different  interactions  between  the  droplet  cluster  and  outside 
environment,  the  bubble  could  be  impermeable  or  permeable  to  heat  and  mass  transfer. 
The  quasi-steady  conservation  equations  of  mass,  species,  and  energy  were  solved  within 
the  bubble.  Results  show  that  saturation  occurs  before  vaporization  can  be  complete  in 
small  bubbles.  Ambient  pressure  and  temperature  produce  large  effects  on  droplet 
interactions.  An  increase  of  either  of  them  could  cause  a  shift  from  incomplete  to 
complete  vaporization. 

Bellan  and  Cuffel  (1983)  developed  a  theory  of  dense  spray  evaporation  based  on 
the  scheme  of  single  droplet  vaporization  in  a  finite  sphere  and  the  globe  conservation 
equations  for  the  two-phase  mixture.  Just  like  Tishkoff  (1979),  the  monosized  droplets 
were  assumed  to  distribute  uniformly  inside  the  droplet  cluster.  Each  droplet  was 
surrounded  in  a  finite  space,  called  “sphere  of  influence”,  the  radius  of  which  is  the  half- 
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distance  between  the  centers  of  two  neighboring  droplets.  The  interstitial  regions  among 
the  spheres  of  influence  possessed  uniform  time-dependent  variables,  which  were 
obtained  by  solving  the  globe  conservation  equations.  A  method  for  solving  single 
droplet  vaporization  was  applied  within  the  sphere  of  influences.  The  variables  in  the 
interstitial  regions  serve  as  the  outside  boundary  conditions.  Since  the  variables  are 
uniform  within  the  entire  interstitial  regions  and  each  droplet  evaporates  in  a  same 
manner,  the  globe  conservation  equations  for  the  droplet  cluster  could  be  simplified  to 
include  only  the  single  droplet,  its  sphere  of  influence,  and  its  related  interstitial  region. 
The  only  difference  between  this  model  and  the  model  proposed  by  Tishkoff  (1979)  is  the 
treatment  of  the  interstitial  region,  which  is  considered  as  a  separated  space  in  this  model, 
but  incorporated  into  the  bubble  in  the  “droplet-in-bubble”  model.  Results  obtained  with 
an  adiabatic  n-decane  droplet  cluster  evaporating  in  the  air  indicated  that  there  are 
distinct  characteristics  that  could  not  be  described  adequately  in  the  single  and  dilute 
droplet  vaporization  theories.  Bellan  and  Harstad  (1987)  further  extended  this  dense 
droplet  vaporization  model  to  include  convection  effect.  In  the  new  method,  a  droplet  still 
evaporated  in  a  quiescent  sphere  of  influence,  while  the  convection  effect  on  the 
vaporization  rate  was  accounted  for  using  an  empirical  correlation. 

All  these  early  investigations  are  related  to  low-pressure  droplet  cluster  behavior. 
Jiang  and  Chiang  (1994)  studied  the  vaporization  of  a  dense  droplet  cluster  in  both  sub- 
and  super-critical  environments.  Droplet  interactions  were  accounted  for  using  the 
“sphere  of  influence”  concept  (Bellan  and  Cuffel  1983).  The  transient  conservation 
equations  for  both  gas  and  liquid  phases  within  the  sphere  of  influence  were  solved  with 
consideration  of  high-pressure  effects,  including  thermodynamic  nonideality  and 
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transport  anomaly,  which  followed  the  theoretical  formulations  of  Hsieh  et  al.  (1991)  for 
the  single  droplet  case.  The  boundary  conditions  at  the  edge  of  the  sphere  were 
determined  by  solving  the  globe  conservation  equations  for  the  interstitial  regions,  where 
space-uniform  variables  were  assumed.  An  adiabatic  and  impermeable  cluster  of  n- 
pentane  droplets  vaporization  in  gaseous  nitrogen  was  studied  over  a  pressure  range  of  1- 
60  bar.  The  predictions  for  dilute  sprays  are  in  qualitative  agreement  with  those  of  single 
droplet  vaporization.  Droplet  lifetimes  are  significantly  prolonged  within  a  dense  spray, 
and  complete  vaporization  can  not  be  achieved  once  a  threshold  value  of  the  sphere  of 
influence  is  reached.  Droplet  interaction  effects  decrease  monotonically  with  the  increase 
of  ambient  pressure. 

The  same  authors  (Jiang  and  Chiang  1996)  improved  their  earlier  method  to 
include  the  interactions  between  the  droplet  cluster  and  the  outside  high-temperature 
environment.  The  sphere  of  influence  was  allowed  to  move  (expansion  or  contraction)  as 
a  result  of  the  constant-pressure  condition.  The  variables  in  the  interstitial  regions  were 
still  assumed  to  be  uniform  in  space,  but  varied  temporally  through  interactions  with  both 
the  sphere  of  influence  and  the  outside  gaseous  environment.  The  outside  gaseous  field 
was  assumed  quiescent  and  spherically  symmetric,  where  the  variable  distributions  were 
determined  by  solving  the  one-dimensional  conservation  equations  of  mass,  energy,  and 
species.  A  dense  cluster  of  n-pentane  droplets  vaporization  in  nitrogen  environments  was 
numerically  investigated  at  two  pressures  of  5  and  60  bar.  It  was  concluded  that,  at  a 
supercritical  condition  (e.g.,  60  bar),  the  surrounding  temperature  outside  the  droplet 
cluster  has  only  limited  effects  on  the  heat-up  of  the  droplet  cluster.  This  is  an  important 
conclusion,  because,  in  this  study,  the  assumption  was  that  the  interstitial  regions  over  the 
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entire  droplet  cluster  were  simultaneously  heated  up  by  the  outside  high-temperature  gas 
environment.  In  practical  situations,  thermal  penetration  should  take  a  longer  time  for  the 
droplets  in  the  inner  region  of  the  cluster  to  feel  the  influence.  Therefore,  heat  transfer 
with  the  outside  environment  should  not  affect  the  droplet  in  the  inner  region  of  the 
droplet  cluster. 

Harstad  and  Bellan  (1998)  developed  a  numerical  model  to  investigate  the  droplet 
behavior  in  a  cluster,  which  is  based  on  the  “sphere  of  influence”  concept  (Bellan  and 
Cuffel  1983)  coupled  with  the  high-pressure  isolated  droplet  formulations  (Harstad  and 
Bellan  1998).  Heat  and  mass  transfer  from  the  outside  environment  to  the  cluster  are 
accounted  for  using  an  empirical  correlation  of  the  Nusselt  number.  A  cluster  of  liquid 
oxygen  (LOX)  droplets  vaporization  in  hydrogen  environments  at  supercritical  conditions 
was  studied,  with  pressure  ranging  from  10  to  80  MPa.  Because  the  value  of  the  Nusselt 
number  was  chosen  arbitrarily,  sensitivity  studies  were  conducted.  Results  indicated  that 
droplet  vaporization  in  the  droplet  cluster  is  insensitive  to  the  value  of  Nusselt  number 
over  three  orders  of  magnitude,  which  is  consistent  with  the  conclusion  of  Jiang  and 
Chiang  (1996)  that  outside  heat  transfer  to  the  droplet  cluster  only  presents  minor  effect. 
1.1.2  Droplet  Vaporization  in  Convective  Environments 

In  practical  combustion  devices,  droplets  are  generally  moving  relative  to  the 
gaseous  environment.  Therefore,  the  study  of  droplet  vaporization  in  high-pressure 
convective  environments  is  a  very  important  research  issue.  Under  convective  conditions, 
a  droplet  may  undergo  severe  deformation  and  even  break-up.  The  droplet  geometry 
becomes  an  unknown  shape,  which  must  be  solved  as  part  of  the  solution.  In  order  to 
solve  this  problem,  a  numerical  algorithm  must  have  the  ability  to  track  an  arbitrarily 
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shaped  moving  interface  or  discontinuity.  For  moving  boundary  problems,  a  few  methods 
have  been  developed  to  track  or  capture  the  moving  interface,  which  can  be  divided  into 
two  categories:  the  front-capturing,  and  front-tracking  approaches.  In  the  front-capturing 
method,  the  entire  computational  domain  is  discretized  using  a  finite  volume  method,  and 
artificial  viscosity  is  possibly  added  to  avoid  oscillations.  The  enthalpy  method  (Voller 
and  Prakash  1987)  used  in  the  material  processing  area  belongs  to  this  category.  The 
other  algorithms  are  the  marker-in-cell  (Harlow  and  Welch  1965),  volume-of-fluid 
(VOF)  (Hirt  and  Nichols  1981),  and  level-set  (Sussman,  Smereka,  and  Osher  1994) 
methods.  In  all  these  treatments,  stationary  Eulerian  grids  are  utilized  and  the  interface 
position  and  shape  are  reconstructed  after  each  computational  step  by  using  markers, 
fractional  volume,  or  a  level-set  function,  which  are  advected  with  the  local  flow 
velocities.  The  disadvantage  of  these  methods  is  that  while  merger  and  break-up  of  the 
interface  can  be  well  handled,  they  can  not  be  treated  with  precision  because  the  interface 
positions  are  not  known  explicitly.  In  the  front-tracking  methods,  the  discontinuous 
interface  is  represented  by  a  set  of  connected  points,  which  are  tracked  explicitly  during 
the  entire  computational  time.  The  other  numerical  grids  can  be  either  stationary,  such  as 
in  the  cut-cell  method  (Udaykumar  1994)  and  the  immersed  boundary  technique 
(Tryggvason  and  Unverdi  1990,  Unverdi  and  Tryggvason  1992),  or  in  motion,  such  as  in 
the  arbitrary  Lagrangian-Eulerian  (ALE)  method  (Hirt,  Amsden,  and  Cook  1974)  and 
arbitrary  moving  mesh  method  (Welch  1993,  Welch  1995).  In  the  immersed  boundary 
method,  the  interface  is  not  kept  sharp  but  rather  given  a  finite  thickness  at  the  order  of 
the  mesh  size.  Thus,  this  method  is  not  suitable  for  solving  the  high-pressure  phase 
change  problems.  The  cut-cell  method  is  initially  developed  to  treat  solid-liquid 
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interfaces  in  the  solidification  problems.  For  the  fluid-fluid  interfaces,  however,  the 
method  is  proved  to  be  sensitive  to  noise  and  is  difficult  to  stabilize.  Therefore,  the 
current  methods  that  are  appropriate  for  solving  high-pressure  phase  change  problems  are 
the  arbitrary  Eulerian-Lagrangian  (ALE)  method,  which  was  applied  for  droplet 
deformation  and  vaporization  by  Deng  and  Jeng  (1992),  and  the  arbitrary  moving  mesh 
method.  In  both  methods,  there  are  explicit  discontinuous  interfaces  where  the  boundary 
conditions  regarding  phase  change  can  be  applied. 

1. 1.2.1  Single  Droplet  Vaporization  in  Convective  Environments 

Welch  (1993,  1995)  studied  a  single-component  liquid  droplet  vaporization  in  its 
own  vapor.  This  is  not  a  droplet  vaporization  case  in  spray  combustors,  but  there  are 
physical  similarities  between  them,  such  as  the  treatment  of  the  phase-changing  interface 
and  the  interfacial  thermodynamics.  In  this  work,  the  interface  was  tracked  explicitly  by 
the  computational  nodes  representing  the  liquid  and  vapor  interfaces.  Two  nodes  located 
in  the  same  spatial  position  are  needed  to  represent  both  phases.  The  interface  motion  is 
found  by  solving  the  interfacial  boundary  conditions  and  the  bulk  flow  equations.  The 
movements  of  the  computational  nodes  in  the  bulk  regions  were  designed  to  limit  the 
amount  of  distortion  of  the  mesh  near  the  interface  during  deformation  process.  Their 
positions  at  the  new  time  level  were  determined  by  interpolation  between  the  movements 
of  the  outside  boundary  and  the  interface  at  the  previous  time  level.  Triangular  meshes 
were  utilized,  which  is  more  endurable  in  terms  of  the  mesh  distortion.  Despite  all  these 
efforts,  the  fluid  flows  are  still  restricted  to  those  in  which  the  interfaces  do  not  undergo 
significant  distortions  or  change  topological  features. 
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Deng  and  Jeng  (1992)  studied  n-heptane  droplet  vaporization  in  nitrogen  gas  at  40 
atm,  with  the  inflow  gas  temperature  and  velocity  ranging  form  500  to  1000  K,  2  m/s  to 
10  m/s,  respectively.  Special  attention  was  devoted  to  droplet  deformation  and  breakup 
using  the  Arbitrary  Lagrangian-Eulerian  (ALE)  method.  In  this  algorithm,  the 
computational  grids  are  first  moving  at  local  velocities.  Then  a  rezone  step,  during  which 
a  new  set  of  grids  is  generated  and  all  the  variables  are  updated,  is  utilized  to  avoid  grid 
distortion.  It  is  very  expensive  to  generate  a  new  set  of  grids  and  update  variables  at  each 
time  step.  Despite  of  all  these  efforts,  the  scheme  still  could  not  handle  flow  situations 
after  droplet  break-up.  Computation  was  terminated  once  the  break-up  occurred.  Results 
indicated  that  the  interfacial  momentum  exchange  between  gas  and  liquid  not  only  moves 
the  droplet  downstream,  but  also  deforms  the  droplet.  Two  types  of  droplet  breakup,  bag- 
and  stripping-types,  were  found,  which  are  sensitive  to  the  gas  phase  velocity  but 
insensitive  to  the  gas  temperature. 

Because  of  the  lack  of  reliable  methods  to  track  the  phase  changing  interface  with 
an  unknown  shape,  the  problem  of  forced-convection  droplet  vaporization  is  very 
difficult  to  handle.  In  order  to  avoid  this  difficulty,  one  common  treatment  is  to  assume 
the  droplet  remains  a  spherical  shape  during  the  entire  vaporization  process,  which  is 
valid  only  with  the  presence  of  extremely  strong  surface  tension  at  low  ambient 
pressures.  A  few  works  have  been  conducted  for  investigation  of  the  vaporization  of 
hydrocarbon  fuels  in  convective  environments. 

Haywood  et  al.  (1989)  made  a  detailed  examination  of  the  evaporation  of  an  n- 
heptane  droplet  in  air  at  800  K  and  1  atm  with  an  initial  Reynolds  number  of  100.  The 
effects  of  variable  thermophysical  properties,  liquid  phase  motion  and  heating,  regression 
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of  droplet  size,  and  transient  variations  of  droplet  velocity  were  considered  in  the 
numerical  analysis.  The  droplet  was  assumed  to  maintain  a  spherical  shape  without 
deformation.  Relationships  were  established  for  the  transient  variations  of  heat  transfer 
rate  (Nusselt  number),  mass  transfer  rate  (Sherwood  number),  and  drag  coefficient. 

The  similar  problem  was  also  studied  systematically  by  Chiang  et  al.  (1992), 
where  variable  properties,  regression  of  droplet  size,  transient  droplet  heating,  and 
deceleration  of  the  flow  due  to  the  drag  action  on  the  droplet  were  included.  A  spherical 
droplet  was  also  assumed.  Different  fuel  types,  such  as  n-octane,  n-decane,  n-hexane, 
were  considered,  with  the  initial  Reynolds  numbers  ranging  from  50  to  150.  Other 
parametric  studies  include  initial  droplet  temperature  and  ambient  temperature. 
Correlations  regarding  transient  heat  and  mass  transfer  rates  and  droplet  drag  coefficient 
were  established,  which  are  different  from  those  of  Haywood  et  al.  (1989)  due  to  strong 
surface  blowing.  The  conclusion  was  that  the  drag  coefficient  might  increase  or  decrease 
over  different  portions  of  the  droplet  lifetime,  while  the  Nusselt  and  Sherwood  numbers 
follow  a  same  gently  decreasing  trend. 

For  a  supercritical  droplet  vaporizing  in  convective  environments,  a  reasonable 
simplification  is  to  assume  the  droplet  reaches  the  critical  mixing  state  immediately  after 
being  injected  into  the  combustor.  This  is  an  approximation  for  most  hydrocarbon  fuels, 
but  it  is  actually  true  for  an  oxygen  droplet  vaporizing  in  hydrogen  environments  at 
supercritical  conditions,  as  demonstrated  by  the  results  regarding  an  isolated  droplet 
vaporizing  in  a  quiescent  environment  (Yang  et  al.  1994,  Lafon  1995).  Once  the  droplet 
reaches  critical  state,  enthalpy  of  vaporization  and  surface  tension  vanish,  and  a 
continuous  flow  field  is  retained.  A  few  research  works  have  been  conducted  in  this  area. 
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Lee  et  al.  (1990)  analyzed  the  dispersion  of  a  preheated  fuel  droplet  that  was 
suddenly  set  into  motion  in  a  gaseous  environment  at  the  same  density.  Temperature  was 
assumed  uniform  in  the  entire  flow  field,  and  density  and  physical  properties  were  treated 
as  constants.  The  Schmidt  number  was  taken  to  be  unit.  The  initial  droplet  injection 
process  was  modeled  by  instantly  imposing  a  harmonic  vortex  sheet  at  the  droplet 
surface,  which  implies  that  the  flow  is  totally  induced  by  the  vorticity.  In  this  anlysis,  the 
transient  axisymmetric  stream  function-vorticity  and  species  equations  were  solved  to 
determine  the  evolution  of  the  vorticity  distribution,  species  mixing,  and  droplet 
deformation.  Two  cases  with  the  Reynolds  numbers  of  50  and  200  were  studied.  The 
study  found  that  a  mushroom-like  droplet  shape  is  formed  due  to  strong  droplet 
deformation,  which  is  caused  by  the  evolution  of  the  vorticity  and  a  formation  of  a  ring 
vortex.  The  conclusion  was  that  fuel  vapor  mixing  is  strongly  produced  by  the  vortex 
ring,  and  the  fuel  concentration  gradient  is  higher  around  the  axisymmetric  axis, 
indicating  an  enhanced  entrainment  of  ambient  gas  and  mixing  at  the  spiral  tip.  The 
effects  are  stronger  at  a  higher  Reynolds  number. 

Daou  and  Rogg  (1998)  investigated  numerically  the  combustion  of  a  supercritical 
fuel  droplet  suddenly  set  into  motion  in  a  high-temperature  oxidizing  environment.  In  this 
work,  thermal  conductivity,  viscosity,  the  constant-pressure  heat  capacity,  and  the 
molecular  weight  of  the  gas  mixture,  as  well  as  the  product  of  density  and  mass 
diffusivity  were  taken  as  constants.  The  Burke-Schumann  flame-sheet  model  was 
adopted,  and  the  Lewis  number  was  assumed  as  unit.  Combustion  was  described  using  an 
irreversible  one-step  reaction  mechanism.  Density  was  further  assumed  to  obey  the  ideal- 
gas  equation  of  state.  At  a  Reynolds  number  of  100,  a  large  recirculating  wake  is  formed 
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at  the  back  of  the  droplet.  The  droplet  experiences  a  substantial  deformation,  especially  a 
stretching  perpendicular  to  the  symmetric  axis,  which  was  attributed  to  the  combined 
effects  of  non-uniform  pressure  distribution  and  vorticity.  The  study  was  focused  on  the 
combustion  time,  which  is  defined  as  the  time  required  to  consume  all  the  fuel  mass.  Due 
to  deformation  and  straining  of  isoscalar  surface,  the  dependence  of  combustion  time  on 
the  magnitude  of  heat  and  mass  diffusion  decreases  rapidly  with  the  increase  of  the 
Reynolds  number.  Once  the  Reynolds  number  exceeds  a  few  hundreds,  the  combustion 
time  becomes  independent  of  the  diffusivity,  and  is  proportional  to  the  convective  time 
multiplied  by  the  square  root  of  density  ratio.  The  study  concluded  that  the  combustion 
time  increases  with  increasing  heat  release  as  a  result  of  gas  expansion  and  the  decrease 
of  the  scalar  gradients.  Droplet  deformation  becomes  weaker  at  a  higher  heat  release. 
Because  of  the  simplicities  regarding  thermophysical  properties,  the  results  of 
supercritical  droplet  combustion  and  fuel  mixing  in  these  two  works  (Lee  et  al.  1990, 
Daou  and  Rogg  1998)  should  be  used  with  caution. 

In  light  of  the  apparent  flaws  with  evaluations  of  thermophysical  properties  in  the 
earlier  works,  Hsiao  (1995)  investigated  systematically  the  problem  of  a  supercritical 
oxygen  droplet  vaporization  in  convective  hydrogen  environments,  where  the  pressure 
effects  on  thermophysical  property  evaluations  have  been  carefully  taken  into  account.  In 
order  to  achieve  high  accuracy  of  density  prediction,  the  Benedict-Webb-Rubin  (BWR) 
equation  of  state  was  utilized  with  an  extended  corresponding  state  (EXCST)  principal. 
The  same  principal  was  also  adopted  for  evaluating  all  the  thermodynamic  properties, 
such  as  enthalpy  and  heat  capacity.  Thermal  conductivity  and  viscosity  were  estimated  by 
the  corresponding  states  method  proposed  by  Ely  and  Hanley  (1981,  1983),  and  the 
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pressure  effect  on  the  binary  mass  diffusivity  was  calculated  by  the  Takahashi  (1974) 
scheme.  Droplet  deformation  and  droplet  lifetimes  at  supercritical  pressures  were  studied 
extensively.  Results  indicate  that  ambient  flow  strongly  influences  droplet  evolution  and 
changes  the  droplet  vaporization  rate  and  dynamic  characteristics.  Correlations  of 
aerodynamic  drag  coefficient  and  droplet  lifetime  were  established  as  functions  of 
ambient  pressure  and  the  Reynolds  number.  Due  to  the  complexity  of  the  BWR  equation 
of  state,  the  pressure  effects  on  the  thermodynamic  relationships  in  the  numerical  scheme 
were  not  considered  properly  in  this  work,  which  could  cause  inconsistency.  In  addition, 
the  thermodynamic  properties  of  a  real  fluid  mixture,  such  as  the  partial  mass  or  partial 
molar  enthalpy,  were  not  taken  into  account 
1.1.2.2  Droplet  Interactions  in  Convective  Environments 

In  convective  environments,  droplet  interactions  present  features  distinct  from 
those  in  a  quiescent  environment.  Because  of  the  complexity  of  droplet  interactions  in 
convective  environments,  numerical  investigations  are  focused  on  the  influences  from  the 
immediate  neighboring  droplets  on  droplet  dynamics,  such  as  drag  coefficient  and  droplet 
collision/separation  phenomenon. 

Patnaik  (1986)  made  an  effort  to  investigate  the  interacting  effects  with  two 
evaporating  droplets  moving  in  tandem  with  respect  to  the  free  stream.  The  complete  set 
of  conservation  equations  of  mass,  energy,  and  species  were  considered  with  full  account 
of  internal  circulation  and  transient  heating  within  the  liquid  phase,  and  variable  density, 
but  all  the  other  thermophysical  properties  were  taken  as  constants.  Droplet  deformations 
were  neglected.  The  two  droplets  were  treated  separately,  and  the  downstream  solution  of 
the  leading  droplet  was  used  as  the  inflow  conditions  of  the  trailing  droplet.  In  this  way. 
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the  influence  of  the  leading  droplet  on  the  trailing  droplet  was  studied,  but  the  influence 
of  the  trailing  droplet  on  the  leading  droplet  was  neglected.  Moreover,  the  transient 
variations  of  the  distance  between  the  two  droplets  could  not  be  investigated. 

Raju  and  Sirignano  (1990)  examined  the  same  problem,  but  allowed  the  two 
droplets  moving  relative  to  each  other.  Studies  were  conducted  over  a  limited  range  of 
initial  Reynolds  numbers,  interdrop  spacings,  and  droplet  radius  ratios.  Results  indicated 
that  when  the  interdrop  spacing  is  sufficiently  large,  both  droplets  behave  as  an  isolated 
droplet.  Decreasing  the  interdrop  spacing,  the  trailing  droplet  is  fully  covered  by  the 
wake  of  the  leading  droplet,  resulting  in  the  reduced  drag  coefficient  and  heat  transfer 
rate  for  the  trailing  droplet.  For  equal  sized  droplets,  collision  occurs  for  all  the  Reynolds 
numbers  and  interdrop  spacings  studied.  Decreasing  the  size  of  the  trailing  droplet  to  a 
certain  limit,  the  two  droplets  never  collide.  The  droplet  radius  ratio  at  this  limit  was 
termed  bifurcation  ratio. 

Chiang  (1990)  extended  the  work  by  Raju  and  Sirignano  (1990)  to  include  the 
variable  thermophysical  properties  due  to  the  conclusion  that  the  drag  coefficient  could 
be  overestimated  by  at  least  20%  for  the  constant-property  case.  Extensive  parametric 
studies  were  conducted  regarding  the  effects  of  initial  Reynolds  numbers,  interdrop 
spacings,  droplet  radius  ratios,  different  fuel  types,  as  well  as  ambient  temperatures  and 
initial  droplet  temperatures.  The  study  concluded  that  the  drag  coefficient  and  Nusselt 
and  Sherwood  numbers  of  the  trailing  droplet  are  significantly  lower  than  those  of  the 
leading  droplet.  The  droplet  spacing  could  increase  or  decrease  with  time,  depending  on 
the  initial  Reynolds  number,  initial  interdrop  spacing,  and  droplet  radius  ratio.  Consistent 
with  the  earlier  results  (Raju  and  Sirignano  1990)  at  a  given  initial  interdrop  spacing. 
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there  exists  a  critical  droplet  radius  ratio  (bifurcation  ratio)  below  which  droplet  collision 
can  not  occur.  Furthermore,  the  critical  ratio  increases  when  the  initial  Reynolds  number 
decreases.  Chiang  (1990)  also  studied  three-droplet-in-tandem  vaporization  in  convective 
environments.  Results  indicated  that  the  general  conclusions  drawn  from  two-droplet 
analysis  could  be  applied  to  the  two  neighboring  droplets  in  the  three-droplet  case. 

All  the  studies  above  are  regarding  droplet  interactions  at  low  pressures.  There  is 
no  investigation  being  conducted  for  the  supercritical  droplet  interactions  in  convective 
environments. 

1.2  Research  Objectives 

The  physical  and  chemical  mechanisms  involved  with  high-pressure  droplet 
vaporization  and  combustion  are  extremely  complicated.  Many  intricate  flow,  transport, 
and  combustion  processes  occur  in  various  parts  of  the  liquid  and  gaseous  phase  regions. 
In  view  of  the  inadequacy  of  existing  theories  and  experimental  techniques,  especially  for 
the  liquid  oxygen  (LOX)  vaporization  in  gaseous  hydrogen  case,  a  numerical 
investigation  of  the  problem  is  conducted  here.  Special  efforts  are  devoted  to  consider  the 
pressure  effects  on  thermophysical  properties.  Thermodynamic  nonideality  and  transport 
anomaly  will  be  treated  properly.  A  robust  all-fluid  numerical  algorithm  is  developed,  in 
which  the  real  fluid  effects  are  treated  consistently  with  full  account  of  pressure  effects. 
For  droplet  vaporization  in  convective  environments,  a  small  Mach-number  region  exists 
inside  the  droplet.  Therefore,  the  preconditioning  numerical  algorithm  with  the  dual  time¬ 
stepping  integration  technique  is  applied  to  solve  the  unsteady  flow  field. 
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The  overall  objective  of  this  work  is  to  answer  some  of  the  fundamental  questions 
regarding  liquid-propellant  droplet  vaporization  and  combustion  in  high-pressure 
environments.  Specific  program  objectives  are: 

1 .  to  develop  a  robust  all-fluid  numerical  algorithm,  which  is  capable  of  treating 
the  fluid  fields  over  the  entire  thermodynamic  regime,  ranging  from  dense 
liquid  to  dilute  gas.  A  unified  treatment  of  the  real  fluids  regarding  both 
thermodynamic  property  evaluation  and  numerical  relationship  derivation  is 
developed  based  on  fundamental  thermodynamic  theories,  which  is  valid  for 
any  equation  of  state.  Full  account  is  also  taken  for  the  pressure  effects; 

2.  to  study  droplet  cluster  behavior  at  both  sub-  and  super-critical  conditions.  In 
view  of  the  complexity  of  the  problem,  the  one-dimensional  simplification  is 
accepted,  and  the  droplet  interactions  are  treated  by  droplet  vaporization  in  a 
“bubble”.  This  represents  the  practical  situation  occurring  in  the  inner  region 
of  a  dense  droplet  cluster; 

3.  to  investigate  droplet  vaporization  in  forced  convective  environments  at 
supercritical  conditions.  A  droplet  will  be  assumed  to  reach  the  critical  mixing 
state  immediately  after  being  injected  into  the  combustor.  The  initial  Reynolds 
number  and  pressure  effects  on  droplet  lifetimes,  droplet  deformations,  and 
droplet  dynamics  will  be  explored  in  detail; 

4.  to  examine  droplet  interactions  in  forced-convection  environments  at 
supercritical  conditions.  Two-droplet-in-tandem  vaporization  in  convective 
environments  is  studied  over  a  wide  range  of  initial  Reynolds  numbers, 
pressures,  and  interdrop  spacings.  The  effects  of  droplet  interactions  on 
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droplet  lifetimes,  droplet  deformations,  and  droplet  dynamics  will  be 
investigated. 

The  current  numerical  research  represents  a  series  of  investigations  into  the 
behavior  of  liquid-propellant  droplets  at  elevated  pressures  and  temperatures.  Results  are 
extremely  helpful  for  understanding  the  fundamental  mechanisms  of  droplet  vaporization 
and  combustion.  This  thesis  is  arranged  in  the  following  manner.  The  theoretical 
formulations  for  an  isolated  droplet  vaporization  in  a  quiescent  environment,  including 
the  conservation  equations  and  thermophysical  property  evaluation  methods,  are 
described  in  Chapter  2.  The  numerical  solution  procedures  are  also  described  in  detail. 
Chapter  3  presents  the  results  regarding  droplet  cluster  behavior.  The  studies  are  focused 
on  the  effects  of  ambient  pressures  and  temperatures.  In  Chapter  4,  a  robust  all-fluid 
numerical  algorithm  with  the  preconditioning  technique  is  developed.  Details  of 
numerical  derivations  are  presented.  The  numerical  algorithm  is  applied  to  study 
supercritical  droplet  vaporization  in  convective  environments  in  Chapter  5,  where  useful 
correlations  of  droplet  lifetime  and  drag  coefficient  are  established.  Chapter  6  discusses 
droplet  interactions  with  two-droplet-in-tandem  vaporization  in  supercritical  convective 
environments.  Conclusions  of  the  present  research  and  recommendations  for  future  works 


are  summarized  in  Chapter  7. 


Chapter  2 


THEORETICAL  FORMULATIONS  OF  SUPERCRITICAL 
DROPLET  VAPORIZATION  IN  QUIESCENT  ENVIRONMENTS 

Based  on  the  extensive  literature  survey,  there  are  clearly  distinct  characteristics 
between  high-pressure  droplet  vaporization  and  its  low-pressure  counterpart.  In  order  to 
model  supercritical  droplet  vaporization  phenomenon  accurately,  the  following  features 
must  be  considered  in  the  numerical  treatment:  1)  Transient  conservation  equations  of 
mass,  energy,  and  species  for  both  the  liquid  and  gaseous  phases  have  to  be  established 
with  the  proper  account  of  the  droplet  surface  regression  rate;  2)  Thermodynamic  non¬ 
ideality,  including  the  pressure  dependent  gaseous  dissolution  in  the  liquid  droplet  and 
the  partial  molar  enthalpy  for  calculating  the  enthalpy  of  vaporization,  must  be  included; 
and  3)  Pressure  effects  on  the  thermophysical  properties  must  be  taken  into  account  with 
proper  treatment  of  transport  anomaly  in  the  transcritical  region.  In  this  chapter, 
theoretical  formulations  of  supercritical  droplet  vaporization  in  quiescent  environments 
will  be  established. 

2.1  Conservation  Equations 

A  droplet  is  assumed  to  evaporate  in  a  quiescent  microgravity  environment,  and 
the  flow  field  is  thus  spherically  symmetric.  The  other  reasonable  simplifications  are  as 
follows: 

1)  Pressure  remains  constant  during  the  entire  vaporization  process  due  to  small  flow 
velocity  compared  with  the  acoustic  wave.  This  assumption  has  been  justified  by  Jia 


and  Gogos  (1993); 
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2)  The  Soret  and  Dufour  effects  are  neglected,  which  is  justified  by  Curtis  and  Farrell 
(1992),  and  Lafon  (1995); 

3)  Thermal  radiation  is  neglected. 

Based  on  these  assumptions,  the  transient  conservation  equations  are  presented  in 
the  following  forms: 

1)  Mass  conservation 

—  f  pdV-v\  p{v  -  w)  •  ndA  ■=  Q  (2.1) 

dt 


2)  Momentum  conservation 

Vp  =  0  (2.2) 

3)  Energy  conservation 

f  k  (’)  M  ^  =  -L,  h.  r,i  P'  ■ 


4)  Species  conservation 


—  f  pY^V+i  pYAv-w)-ndA  =  -[  d-ndA 


(2.4) 


where  /  =  1,  - 1 ,  and  N  is  the  total  species  number.  Therefore,  Eq.  2.4  represents 

N -I  species  equations.  The  physical  variables  p,v,ej,p,Yi  are  density,  velocity, 
specific  total  internal  energy,  pressure,  and  species  concentration  respectively,  vv 
represents  the  moving  velocity  of  the  control  surface.  V„,  A,  are  the  volume  and  surface 

area  of  the  control  volume,  n  is  the  unit  vector  normal  to  the  control  surface.  The 
specific  total  internal  energy  is  expressed  as 

e,  —  e  - 

2 


(2.5) 
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where  eis  the  internal  energy. 

The  diffusive  flux  terms  q^,q.  in  Eqs.  (2.3)  and  (2.4)  are  determined  by  Fourier’s 

and  Pick’s  Laws,  respectively. 

1)  Species  diffusion 

q^=pY,(v,-v)  =  -pD,„yY,  (2.6) 

where  v,.  is  the  velocity  of  species  i.  is  the  effective  mass  diffusivity  of  species  i  in 
the  mixture,  which  is  determined  in  terms  of  binary  mass  diffusivity  D..  (Bird  et  al. 

1960), 

N 

(2.7) 

i^j 

2)  Energy  diffusion 

Energy  diffusion  consists  of  two  parts.  One  is  the  thermal  conduction  term,  and 
another  term  is  related  to  mass  diffusion. 

qe=-^^T-Y^qihi  (2.8) 

i=\ 

where  X  is  thermal  conductivity,  and  /?,•  is  the  partial  mass  enthalpy  of  species  i,  which 

is  defined  later  in  this  chapter. 

2.2  Boundary  Conditions 
2.2.1  Initial  Boundary  Conditions 

Inside  the  liquid  droplet,  initially  a  uniform  temperature  is  prescribed,  and  there  is 
no  fluid  flow  and  no  gaseous  species, 
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(2.9) 

=  0 

^/,0  =  0 

where,  7)  (,,7,,,  are  the  initial  liquid  temperature  and  the  prescribed  liquid  temperature 
respectively.  initial  fuel  mass  fraction  and  the  mass  fraction  of  the 

gaseous  species,  v,  is  the  initial  liquid  velocity. 

The  initial  gaseous  temperature  is  also  assumed  to  be  uniform,  and  there  is  no 
velocity  and  no  fuel  species  in  the  gaseous  phase. 


^«.0 

Y;.o  =  0 
=  1 
v..o=0 


(2.10) 


where  the  subscript  g  refers  to  the  gaseous  phase,  and  all  the  other  nomenclatures  are  the 

same  as  those  in  the  liquid  phase,  as  defined  above. 

At  droplet  surface,  the  conventional  assumption  of  thermodynamic  phase 
equilibrium  holds,  resulting  in  the  equality  of  temperature,  pressure,  and  chemical 
potentials  between  the  two  phases.  Accordingly,  the  equilibrium  mass  fractions  in  both 
sides  can  be  calculated.  The  thermodynamic  phase  equilibrium  boundary  conditions  will 
be  established  in  detail  later. 

2.2.2  Boundary  Conditions  at  Droplet  Center  and  Outside  Field 

Spherical  symmetry  boundary  conditions  remain  at  the  droplet  center. 


vr  =  o 
vy.  =  0 

V  =0 


(2.11) 
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(a) 


Liq  uid 


Vapor 


[p(v  -  R)-nA],. 


Interface  with 
velocity  ^ 


Differential  Control  Volume 
Moving  with  the  Interface 


(b) 


Liquid  Vapor 


velocity  R  M  oving  w  ith  the  Interface 

(C) 


Liquid  Vapor 


velocity  ^  Moving  with  the  Interface 


Fig.  2.1  Boundary  conditions  at  droplet  surface  (a)  mass  conservation; 
(b)  species  conservation;  (c)  energy  conservation 
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Outside  boundary  is  maintained  at  a  distance  of  100  times  of  the  droplet  radius 
from  the  droplet  center.  All  the  variables  there  will  not  be  affected  by  droplet 
vaporization  and  remain  at  their  initial  values. 

2.2.3  Boundary  Conditions  at  Droplet  Interface 

Phase  change  occurs  at  droplet  surface,  where  undergo  rapid  heat  and  mass 
transfer.  Mass,  energy,  and  species  across  the  surface  must  be  conserved,  as  shown  in 
Fig.  2.1.  The  conservation  equations  are 

m  =  =  piv  -  R)-nA\^^^  (2.12) 

m,.  =  [mYi  +  qi  ■  =  [mY.  +  q.  ■  (2. 13) 

-AVr|  (2.14) 

(=1 

where  R  is  the  droplet  regression  rate,  m  the  mass  vaporization  rate,  and  m,-  the 
vaporization  rate  of  species  i. 

Except  for  these  conservation  equations,  thermodynamic  phase  equilibrium  is  also 
assumed  at  the  droplet  surface,  which  means  temperatures,  pressures,  and  chemical 
potentials  at  both  sides  of  the  surface  must  be  equal  (Callen  1985), 

p^=p'  (2.15) 

=Pi'i 

This  general  expression  of  thermodynamic  phase  equilibrium  can  naturally  take 
into  account  the  dissolution  of  gaseous  species  in  the  liquid  droplet,  which  becomes 
extremely  important  with  the  increase  of  pressure.  As  a  specific  example,  the  phase 
diagrams  of  liquid  oxygen  (LOX)  in  gaseous  hydrogen  are  illustrated  in  Fig.  2.2.  In  Fig. 
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2.2,  the  pressure  dependent  solubilities  are  clearly  illustrated,  where  the  critical  mixing 
states  are  also  shown  to  depend  strongly  on  pressure. 


Fig.  2.2  General  phase  diagrams  of  oxygen  and  hydrogen  at  different  pressures 
2.3  Irreversible  Thermodynamic  Phenomenon 

In  this  supercritical  droplet  vaporization  problem,  the  internal  boundary  at  the 
droplet  surface  renders  numerical  treatment  very  time-consuming.  In  fact,  the  boundary 
conditions  established  above  are  still  not  sufficient  for  solving  this  problem,  and  one 
more  boundary  condition  is  required.  In  the  early  works  reviewed  in  Chapter  1,  two 
methods  are  employed  to  handle  this  problem,  which  are  directly  related  to  the 
characteristics  of  spherical  symmetry.  One  is  to  assume  the  velocity  within  the  liquid 
droplet  is  zero.  This  is  a  reasonable  assumption  because  the  liquid  velocity  is  indeed  very 
small  without  internal  circulation.  Another  one  is  to  relate  the  interfacial  regression  rate 
to  the  total  mass  conservation  in  the  liquid  phase. 
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In  this  study,  a  more  efficient  method  is  utilized,  which  is  initially  proposed  by 
Lafon  (1995)  based  on  irreversible  thermodynamic  theory.  According  to  Onsager’s 
irreversible  theorem  (Callen  1985),  all  the  irreversible  fluxes,  including  heat  and  mass 
fluxes,  can  be  related  to  some  generalized  forces  called  affinities,  which,  in  simple 
thermodynamic  systems,  include  the  differences  of  temperature,  pressure,  and  chemical 
potential  of  each  species.  In  this  theorem,  the  dependence  of  a  flux  on  the  affinities  is 
nonlinear.  However,  the  high  order  terms  are  always  found  to  be  negligible.  In  fact,  as 
Fourier’s  law  indicates,  heat  flux  is  only  strongly  affected  by  temperature  gradient,  which 
is  the  driving  force  for  heat  transfer.  Thermodynamics  indicates  that  the  driving  force  for 
species  transfer  is  chemical  potential.  Moreover,  thermodynamic  phase  equilibrium 
assumption  requires  extremely  small  differences  of  temperature  and  pressure  at  the 
droplet  surface.  Therefore,  it  is  physically  reasonable  to  neglect  the  effects  of  temperature 
and  pressure  on  mass  transfer,  which  leads  to  the  following  linear  correlation; 

(2.16) 

This  relation  is  extremely  useful  for  a  simple  interfacial  boundary  condition 
treatment.  The  only  problem  is  how  to  determine  the  eonstant  . .  In  general,  there  is 

no  specific  formula  available  to  determine  this  constant.  However,  because 
thermodynamic  phase  equilibrium  is  expected  at  the  droplet  surface,  this  constant  can  be 
simply  chosen  as  a  very  large  number  to  satisfy  our  requirement  of  small  differences  of 
chemical  potentials.  The  accuracy  of  this  method  has  been  demonstrated  by  Lafon  (1995) 
using  results  from  an  isolated  droplet  vaporization  in  quiescent  environments. 
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2.4  Equations  of  State 

In  order  to  solve  the  thermodynamic  phase  equilibrium  problem,  as  well  as  to 
evaluate  thermophysical  properties,  an  equation  of  state  capable  of  handling  real  fluid 
behavior  and  high-pressure  effect  is  required.  The  most  common  equation  of  state 
utilized  for  calculating  high-pressure  phase  equilibrium  is  the  cubic  equations  of  state, 
which  includes  Redlich-Kwong,  Peng-Robinson,  and  the  modified  Soave-Redlich- 
Kwong  (SRK)  equations  of  state.  The  modified  SRK  equation  of  state,  which  is  capable 
of  handling  the  quantum  gas  behavior  of  hydrogen,  is  adopted  here  for  the  treatment  of 
thermodynamic  phase  equilibrium,  calculation  of  thermodynamic  properties,  and 
derivation  of  numerical  relationships.  The  consistent  treatment  of  thermodynamic 
properties  using  one  equation  of  state  produces  an  efficient  numerical  algorithm.  A  more 
complicated  Benedict-Webb-Rubin  (BWR)  equation  of  state,  combined  with  the 
extended  corresponding  states  theory,  is  utilized  for  estimating  transport  properties,  as 
suggested  by  Ely  and  Hanly  (1981,  1983).  Both  the  modified  SRK  and  the  BWR 
equations  of  state  are  capable  of  representing  liquid  phase  behavior.  The  BWR  equation 
of  state  can  be  applied  over  broader  temperature  and  pressure  ranges,  but  the  SRK 
equation  of  state  is  easier  to  apply,  especially  for  deriving  the  thermodynamic  differential 
relationships. 

2.4.1  Modified  Soave-Redlich-Kwong  (SRK)  Equation  of  State 

The  modified  SRK  equation  of  state  is  considered  to  be  both  simple  and  fairly 
accurate,  which  takes  the  following  form  (Graboski  and  Daubert  1978,  1979): 

aa 

^~{W-bp)  ~W  {W  +  bp) 


(2.17) 
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where  is  the  universal  gas  constant.  The  parameters  ‘a’  and  ‘b’  account  for  the  effects 

of  the  attractive  and  repulsive  forces  between  molecules,  respectively,  ‘a  ’  is  the  third 
parameter,  which  is  a  function  of  temperature  and  acentric  factor.  For  mixtures,  they  are 
calculated  from  the  following  mixing  rules: 

N  N 

ecu  —  ^  '^^XjX jCCjjUjj  (2.18) 

/=l7=l 

N 

b=Y.Xibi  (2.19) 

(=1 

The  cross  parameter  a.ja..  in  Eq.  (2.18)  is  given  by 
N  N 

~  ^  ^ ^ J  ~  ^ij  ^  (2.20) 

/=l7=l 

where  x.  is  the  mole  fraction  of  species  i.  k.j  is  the  binary  interaction  coefficient,  which 

can  be  calculated  from  a  known  data  set  (Graboski  and  Daubert  1978).  The  constants 
ai,b.  can  be  determined  from  the  following  universal  relationships: 

2  2 

a- =0.42747  (2.21) 

Pci 

bj  =0m664^^  (2.22) 

Pci 

The  third  parameter  for  species  i  is  given  as 

a,=[l  +  5,,(l-Vr)r  (2.23) 

where  T^,,  p^,,  are  the  critical  temperature  and  pressure  for  species  i.  is  the  reduced 
temperature  for  species  i,  and  S.  is  another  variable.  These  parameters  are  given  by  the 


following  formulas. 


5.  =0.48508 +  1.55 17(y,.  -0.15613ffl' 


(2.25) 


Equations  (2.23)  and  (2.25)  are  not  capable  of  accurately  correlating  hydrogen 
systems  due  to  its  quantum  gas  behavior.  A  modified  expression  of  a  for  hydrogen  is 
further  established  (Gradoski  and  Daubert  1979), 

=  1 .202  exp(-0.30228r )  (2.26) 

This  correlation  is  expected  to  be  accurate  for  hydrogen  at  reduced  temperatures 
greater  than  about  2.5  (the  critical  temperature  of  hydrogen  is  33.2K).  No  mixture  of 
technical  interests  lies  below  this  temperature.  When  Eq.  (2.26)  is  used,  the  binary 
interaction  coefficients  K.j  involving  hydrogen  species  should  be  set  to  0. 

2.4.2  Benedict-Webb-Rubin  (BWR)  Equation  of  State 

The  BWR  equation  of  state  of  methane  is  presented  here,  which  is  proposed  by 
Jacobsen  and  Stewart  (1973)  based  on  very  broad  ranges  of  experimental  data  to  ensure 
the  proper  extrapolation  to  low-temperature  and  high-pressure  regions.  Its  expression  is 

p  =  ^A^{T)p^  +  (2.27) 

n=l  n=10 

where  the  15  coefficients  are  given  in  Table  2.1,  which  are  functions  of  temperature. 

The  32  semi-empirical  Constants  of  Ns  are  list  in  Table  2.2.  The  units  for  pressure, 
density,  and  temperature  are  given  in  bar,  mol/liter,  and  K.  The  stain-rate  y  is  0.0096 
(mol/liter)'^. 


Table  2.1  Coefficients  of  An 


A|=RuT 

A9=  Ni9/T^ 

A2=NiT+N2T'^^+N3+N4/T+N5/T^ 

Aio=N2o/T'+  N2i/T^ 

A3=N6T+  N7+  N8/T+  N9/T^ 

All=  N22/T^+  N23/T^ 

A4=  NioT+  N]i+  N12/T 

Ai2=N24/T'+  N25/T' 

A5=N,3 

Ai3=  N26A^"+  N27/T^ 

A6=  N,4Ar+  n,5/t' 

Ai4=  N28/T^+  N29/T^ 

A7=  NieAT 

A8=  N17/T+  N,8/T^ 

Ai5=  N3(/r"+  N32/T' 

Table  2.2  Semi-empirical  Constants  Based  on  Methane 


Ni=-1.184347314485E-2 

N2=7.540377272657E-1 

N3=-1.225769717554E+1 

N4=6.260681393432E+2 

N5=-3.490654409121E+4 

N6=5.301046385532E-4 

N7=-2.875764479978E-1 

N8=5.011947936427E+1 

N9=-2.821562800903E+4 

Nio=-2.064957753744E+5 

Ni  1=1. 28595 1844828E-2 

Ni2=-1.106266656726E+0 

Ni3=3.060813353408E-4 

Ni4=-3.174982181302E-3 

Ni5=5.191608004779E+0 

Ni6=-3.074944210271E-4 


Ni7=1.071143181503E-5 

Ni8=-9.290851745353E-3 

Ni9=1.610140169312E-4 

N2o=3.469830970789E+4 

N2i=-1.370878559048E+6 

N22=1.790105676252E+2 

N23=1.615880743238E+6 

N24=6.265306650288E-1 

N25=1.820173769533E+1 

N26=1.449888505811E-3 

N27=-3.159999123798E+1 

N28=-5.290335668451E-6 

N29=1. 694350244 152E-3 

N3o=8.612049038886E-9 

N3i=-2.598235689063E-6 

N32=3.153374374912E-5 
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2.5  Thermodynamic  Properties 

In  order  to  solve  the  conservation  equations,  additional  thermodynamic 
relationships  relating  thermodynamic  properties  to  temperature  and,  in  high-pressure 
cases,  pressure,  are  required.  The  thermodynamic  properties  to  be  evaluated  are  density, 
internal  energy,  enthalpy,  entropy,  partial  molar  enthalpy  and  chemical  potential  of  each 
species.  These  properties  can  be  derived  directly  from  fundamental  thermodynamic 
relationships. 


e{T,p)  =  eQ{T)  +  f^^ 


P  P 


dp 


Jr 


(2.28) 


h(T,p)  =  ho{T)+‘ 


PO 


l^T_{dp^ 

P  p'^ 


dT 


yui  Jp 


dp 


(2.29) 


s{T,p)  =  sq{T,p^)-\P 


P^ 


^  dp  ^ 


dp 


Jr 


(2.30) 


Cp(T,p)^Cyo(T)-j^ 

^P( 


P 


dp  +  ■ 


r  On  P  ^ 


p 


(  ^ 

ydpjT 


(2.31) 


where  the  e,h,s,Cp  Cy  are  internal  energy,  enthalpy,  entropy,  constant-pressure  heat 


capacity,  and  constant-volume  heat  capacity,  respectively.  The  subscript  0  refers  to  an 
ideal  state  at  a  low  pressure.  All  the  partial  derivatives  in  these  relations  are  calculated 
based  on  the  Soave-Redlich-Kwong  equation  of  state,  and  they  are 


dT 


PK 


1 


JPj 


,.  {M^-bp)  M 


w 


dT 


(a  a) 


P 


?,F/  (M  w ^P) 


(2.32) 


M^R^T  gq  p(2My^+bp) 
{M^-bpf  {M^+bpf 


dp 


JT 


(2.33) 
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^  dp  ^ 


- [M^  +  p{bi-b)] 


M^.{M^~bpf 


2pJ^Xjaija,j 


■  + 


9 

aap  bf 


^  Wi  (M  W  ^p)  M  yy  -  (M  yy  +  bp^' 


(2.34) 


where  the  derivative  — {aa)  is  given  in  Appendix  B. 
oT 


As  being  demonstrated  in  the  early  works  (Manrique  and  Borman  1969,  Hsieh  et 
al.  1991),  the  properties  of  each  component  in  a  non-ideal  mixture  have  to  be  defined  by 
the  partial  molar  or  partial  mass  properties  in  order  to  find  faithful  results.  According  to 
thermodynamics,  any  thermodynamic  property  in  a  mixture  is  a  function  of 
temperature,  pressure,  and  the  masses  of  the  constituent  species, 

m(/)  =  m0{T,  p,mi)  (2.35) 


Therefore,  the  partial  mass  property  is  defined  as 
dm(p 


<I>i 


V  JT,p,m jjti 


(2.36) 


where  /,  y  =  1,  •  •  • ,  ,  and  ^  refers  to  any  proper  thermodynamic  property  per  unit  mass 

of  a  mixture,  such  as  enthalpy  and  internal  energy.  Figure  2.3  presents  the  variations  of 
the  enthalpy  of  vaporization  of  oxygen  at  thermodynamic  phase  equilibrium  states,  where 
the  enthalpy  of  vaporization  is  calculated  based  on  the  partial  mass  enthalpy  and  shows 
strong  dependence  on  both  temperature  and  pressure. 

In  fluid  mechanics,  density  is  more  convenient  to  use  instead  of  mass.  Therefore, 
it  is  better  to  define  the  mixture  thermodynamic  properties  per  unit  volume. 


p(p^p^{T,p_) 


(2.37) 
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Temperature  (K) 


Figure  2.3  Variations  of  the  enthalpy  of  vaporization  of  oxygen  at  thermodynamic  phase 
equilibrium  states. 

This  leads  naturally  to  the  definition  of  the  partial  density  property  ^ , 


<Pi  = 


(2.38) 


According  to  this  definition,  the  partial  density  internal  energy,  enthalpy,  and 
entropy  are  expressed  as 


e;  = 


^dpe'^ 


(2.39) 


K  = 


^dph'^ 

V  JT.PJ^, 


(2.40) 


^  dps^ 
V  ^P‘  J 


P’Pj^i 


(2.41) 
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where  calculations  of  and  s.  require  the  thermodynamic  relations  presented  in  Eqs. 

(2.28-2.30),  and  the  SRK  equation  of  state.  In  addition,  a  relationship  between  the  partial 
density  property  and  the  partial  mass  property  is  derived  as 


(2.42) 


Utilizing  fundamental  thermodynamic  theories,  including  Euler  equation  and 
Gibbs-Duhem  relation,  the  following  relationship  about  chemical  potential  of  each 
species,  which  is  required  for  the  treatment  of  thermodynamic  phase  equilibrium,  can  be 
derived: 


^■=f.=ei-Tsi  (2.43) 

The  details  for  calculating  thermodynamic  properties  are  presented  in  Appendix 
A.  As  a  specific  example,  the  partial  density  internal  energy  and  enthalpy  of  species  i ,  e) 


ei  can  be  expressed  as 
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Substituting  Eq.  (2.42)  into  Eq.  (2.48b),  and  taking  use  of  the  fundamental 
expressions  of  enthalpy,  which  are  available  in  any  thermodynamics  textbook,  the 
following  relationship  regarding  the  partial  mass  enthalpy,  hj ,  can  be  derived: 
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^  dp  ^ 


Aj  Cj  + 


T 

yui 

^Pj 

P] 

'dp] 

ydpjTY 


dp  ^ 
dpi 


JT  ,/? 


(2.49) 


2.6  Transport  Properties 

Transport  properties  to  be  evaluated  include  viscosity,  thermal  conductivity  and 
binary  mass  diffusivity.  Accurate  estimations  of  these  properties  are  very  important  for 
high-pressure  droplet  vaporization  computation,  since  they  will  determine  the  fluid  field, 
as  well  as  heat  and  mass  transfer  rates.  As  proved  in  the  literature,  pressure  effects  have 
to  be  taken  into  account  to  estimate  these  properties. 

2.6.1  Corresponding  States  Theories 

The  law  of  corresponding  state,  which  was  originally  proposed  by  van  der  Waals 
in  1873,  expresses  that  the  equilibrium  properties  of  different  fluids  can  be  related  to  their 
critical  properties  in  a  universal  way.  For  example,  if  pressure,  temperature,  and  volume 
are  related  to  their  critical  properties,  the  PVT  function  relating  the  reduced  pressure, 
temperature,  and  volume  becomes  identical  for  all  substances.  The  reduced  property  is 
defined  as 

p,  =  pi  Pc,T,=T/T^,V,=V /Vc  (2.50) 

The  corresponding  states  theory  also  holds  for  other  properties,  including 
viscosity,  and  thermal  conductivity.  Via  the  corresponding  states  argument,  a  property  of 
any  fluid  can  be  estimated  by  relating  to  its  counterpart  of  a  reference  substance,  whose 
property  can  be  easily  determined  (Ely  and  Hanly  1981). 

nx(T'’P)  =  Vo(To,Po)Fjj  (2.51) 


where  the  subscript  x  refers  to  the  fluid  of  interest,  0  to  the  reference  fluid. 
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- 


Mw^ 

Mwq 


sl/2 


1/2, -2/3 


Jx,0  %,0 


(2.52) 


and  where  Mw  is  the  molecular  weight.  To  and  po  of  the  reference  fluid  are  calculated  by 

(2.53) 


T(,=T/f,^0 
Po  ~  P^x,0 


In  general  corresponding  states  theory,  the  parameters  f^  Q  and  hj^  Q  is  the  ratio 

of  the  critical  temperature  and  volume.  However,  the  range  of  its  applicability  can  be 
broadened  considerably  by  introducing  the  extended  corresponding  states  (EXCST) 
model,  where  the  equivalent  parameters  of  f^  Q  and  h^  Q  become 

fx,0  ~  ^c,x  ^^c,o')^^r,x’^r,x>^x')  (2.54) 

^A',0  ~  (^C,A  (^C,0^(2r,A’^r,A’ )  (2.55) 

In  Eq.  (2.54)  and  (2.55),  6  and  (f)  are  the  so-called  shape  factors,  which  are  functions  of 
Pitzer’s  acentric  factor  co^  and  the  reduced  temperature  and  volume.  These  functions  can 


be  calculated  by  the  following  correlations: 

^(2r,A'^r,A’^A)“  ^  (^A  “  ^o)-^(2r,A>^r,A) 

^^r,x'^r,x'^x^~  (®A  ~  ^o)^(2"r,A’^r,A)J^c,0  ^ ^c,> 

where  Z  is  the  compressibility  factor. 


(2.56) 

(2.57) 


^(7’r,A>^r,A)=«i +^iin7’; +(ci + Ji/t;)(v;  -0.5) 


G{Tr,, yr,x)  =  Cl2  (^A^  +  ^2 )  +  ^2  (^a""  +  ^2  ) In 


(2.58) 


(2.59) 


and  where 


=  min[2,  max(7).  ^(.,0.5)] 


(2.60) 


=  min[2,max(y^  ;^,0.5)] 

The  coefficients  in  Eqs.  (2.58)  and  (2.59)  are  list  in  Table  2.3. 
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(2.61) 


Table  2.3  Coefficients  for  Shape  Factor  Correlations 


Coefficients  in  Eq.  (2.58) 

Coefficients  in  Eq.  (2.59) 

ai=0.090569 

ai=0.394901 

Bi=-0.862762 

bi=-1.023545 

ci=0.316636 

ci=-0.932813 

Di=-0.465684 

di=-0.754639 

The  extended  corresponding  states  theory  will  be  utilized  with  appropriate  mixing 
rules  to  estimate  the  transport  properties  of  mixtures. 

2.6.2  Mixing  Rules  for  Mixtures 

In  order  to  apply  the  corresponding  states  theory  to  a  mixture,  the  variables  of  that 
mixture  have  to  be  calculated  from  the  corresponding  parameters  of  its  components  using 
appropriate  mixing  rules.  The  following  mixing  rules  are  defined  for  calculating  transport 
properties  (Ely  and  Hanly  1981); 


'■  j 


(2.62) 


fxfi  ~ 


I  J 


The  molecular  weight  of  the  mixture  can  be  found  by 


i  j 


(2.63) 


(2.64) 


where  the  indexes,  i  and  j,  represent  the  species  components  in  the  mixture.  The 
following  combining  rules  are  chosen; 

fij,0  = 


h 


(2.66) 
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Mwij  =  IMwiMwj  l{Mwi  +  Mwj)  (2.67) 

where  the  variables,  and  Ijj ,  are  the  binary  interaction  parameters  with  values 


close  to  zero,  which  are  set  to  zero  for  estimating  transport  properties. 

2.6.3  Viscosity  of  Mixtures 

The  extended  corresponding  states  one  fluid  model  (Ely  and  Hanly  1981)  is 
applied  here  for  estimating  the  voscosity  of  a  mixture.  The  basic  idea  is  straightforward. 
First,  it  is  assumed  that  the  viscosity  of  a  single-phase  mixture  is  equated  to  that  of  a 
hypothetical  pure  fluid, 

fl„iAT,p)=MAT,P)  (2-68) 


where  the  subscripts,  mix  and  x,  refer  to  the  mixture  and  the  hypothetical  pure  fluid, 
respectively. 

The  corresponding  states  theory  for  transport  properties  is  then  used  to  evaluate 
the  viscosity  of  the  hypothetical  pure  fluid  with  respect  to  a  given  reference  fluid, 

juAT,P)  =  Mo(To^Po)Ff^  (2.69) 

where  the  subscript,  0,  refers  to  the  reference  fluid. 


p  n1/2 


Mwq 


All, -2/3 

Jx,0  %,0 


(2.70) 


The  variables,  fxfi,hx,0^  Mw^dxt  then  evaluated  by  the  extended 

corresponding  states  (EXCST)  theory  and  the  mixing  rules. 

Methane  is  chosen  as  the  reference  fluid,  because  of  the  existence  of  its  PVT  and 
viscosity  data  over  the  entire  range  of  fluid  states.  The  date  fitted  methane  viscosity 


correlation  is 
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fiQ{pQ,TQ)  -  Pq\Tq)  +  P^q\Tq)Pq  +  ^Pq{Pq,Tq)X  ^  (2.71) 

where  the  first  two  terms  represent  the  dilute  gas  and  first  density  correction, 
respectively,  while  the  last  term  is  a  remainder  which  dominates  in  high-density  region. 
The  factor,  is  a  correction  of  the  possible  non-correspondence  of  viscosity.  The 

details  about  this  correlation  can  be  found  in  the  paper  of  Ely  and  Hanly  (1981). 

2.6.4  Thermal  Conductivity  of  Mixtures 

Estimating  the  thermal  conductivity  of  a  mixture  is  more  complicated,  since  it  is 
affected  by  two  factors,  one  arising  from  the  transfer  of  energy  from  pure  collision  or 
translation  effect,  X ,  and  another  from  the  transfer  of  energy  via  the  internal  degree  of 
freedom,  X  (Ely  and  Hanly  1983).  The  latter  term  is,  in  general,  independent  of  density. 
Therefore,  only  the  collision  or  translation  part  is  evaluated  using  the  extended 
corresponding  states  one  fluid  model,  which  takes  the  same  procedure  as  that  for 
estimating  the  viscosity  of  a  mixture.  First,  the  thermal  conductivity  of  a  mixture  is 
equated  to  that  of  a  hypothetical  pure  fluid  as 

AX,ipJ)  =  r,{p,T)  (2.72) 

The  corresponding  states  theory  is  then  applied  to  calculate  the  thermal 
conductivity  of  the  hypothetical  pure  fluid  via 

4(r,p)  =  ^(ro,Po)^x^x  (2.73) 

where  takes  the  same  form  as  in  Eq.  (2.70),  whose  calculation  requires  the  extended 
corresponding  states  theory  and  the  mixing  rules,  and  is  the  correction  term  for  non- 


comespondence. 
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Finally,  the  thermal  conductivity  of  a  mixture  can  be  expressed  as 

KvAP.T)  =  +  K,uiT) 


In  Eq.  (2.74)  the  thermal  conductivity  of  the  reference  fluid,  methane,  is  evaluated 
via  an  empirical  correlation.  The  second  term  in  Eq.  (2.74),  which  is  related  to  internal 
degree  of  freedom,  can  be  calculated  via  the  modified  Eucken  correlation  with  mixing 
rules  (Ely  and  Hanly  1983),  which  are 


A-Mwj 

pf 


=  1.32(C*,-^^) 


(2.75) 


and, 

K,iAT)  =  2'Z^iXj^  (2.76) 

i  j 


(2.77) 


where  the  indexes,  i  and  j,  refer  to  the  components  in  the  mixture,  ju^  is  the  dilute  gas 


viscosity  of  component  i,  ■  is  the  ideal  gas  heat  capacity  of  component  i,  and  R^^  is 

the  universal  gas  constant. 

2.6.5  Binary  Mass  Diffusivity 

Binary  mass  diffusivity  presents  more  challenge  than  viscosity  and  thermal 
conductivity  because  of  the  following  two  reasons.  1 .  There  are  only  limited  experimental 
data  existent  for  binary  mass  diffusivity  at  high  pressures,  which  results  in  few  estimation 
methods.  2.  There  is  no  satisfactory  liquid  state  theory  available  for  calculating  binary 
mass  diffusivity  in  the  liquid  phase. 
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In  this  work,  the  binary  mass  diffusivity  in  gaseous  phase  at  low  pressure  is 
evaluated  using  empirical  correlation  of  Fuller  et  al.,  which  is  recommended  by  Reid  et 
al.  (1987).  It  takes  the  following  form: 


0.00143r^-^^ 


(2.78) 


where  Dy  is  binary  mass  diffusivity,  cm  /s 


T  is  temperature,  K 
p  is  pressure,  bar 

Mwij  is  the  combined  molecular  weight,  as  expressed  in  Eq.  (2.67) 


and  'Ey  is  found  for  each  component  by  summing  atomic  diffusion  volumes, 
which  is  tabulated  in  Reid  et  al.  (1987,  Table  11-1). 

High-pressure  effect  on  binary  mass  diffusivity  is  evaluated  by  the  method 
proposed  by  Takahashi  (1974),  which  is  based  on  a  simple  corresponding  state  method. 


DjjP 

(Dijpf 


f{Tr,Pr) 


(2.79) 


where  the  superscript  +  indicates  the  low-pressure  values  given  by  Eq.  (2.78).  The 
function  f{Tj.,Pj.)  represents  a  scaling  factor  of  pressure  based  on  the  reduced 
temperature  and  pressure,  which  is  tabulated  by  Takahashi  and  also  shown  in  Reid  et  al. 
(1987,  Fig.  11-3).  In  order  to  calculate  the  reduced  parameters,  the  following  combing 
rules  for  pseudocritical  properties  of  a  mixture  are  used: 

(2.80) 


Pc=^iPc,i+^jPcJ 


(2.81) 
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The  binary  mass  diffusivities  in  liquids  are  evaluated  by  the  method  of  Tyn  and 
Calus,  as  recommended  by  Reid  et  al.  (1987).  For  a  binary  mixture  of  solute  i  in  solvent), 
it  is 


.  V/6 


Dy  =8.93x10 


-8 


V/ 

V  J  J 


/  p  \0.6 


T_ 


(2.82) 


where  V  is  the  molar  volume  at  the  normal  boiling  temperature,  cm  /mol 
T  is  temperature,  K 
jUj  is  the  viscosity  of  solvent,  cP 


Pi  and  Pj  are  parachors  for  the  solute  and  solvent,  and  the  calculation  method  is 


given  in  Reid  et  al.  (1987). 

There  is  no  correlation  currently  available  for  estimating  the  pressure  effect  on 
binary  mass  diffusivity  in  liquid  phase. 

2.7  Numerical  Method 

An  implicit  finite  volume  scheme  is  used  in  this  study  to  solve  the  problem  of 
supercritical  droplet  vaporization  in  a  quiescent  environment.  The  implicit  scheme  is  used 
because  there  is  no  numerical  convergence  restriction  for  time  steps,  which  allows  us  to 
choose  them  purely  from  the  physical  consideration,  leading  to  large  time  steps  and 
reduced  total  computational  time.  The  numerical  methods  are  presented  in  detail  in  this 


section. 
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2.7.1  Grid  Generation 

Because  two  phases  have  to  be  treated  simultaneously,  finite  volumes  will  be 
generated  for  each  phase,  namely  liquid  oxygen  and  gaseous  hydrogen  phases,  separately. 
The  control  surfaces  are  concentric  spheres  with  the  droplet  center  as  the  sphere  center. 
Because  the  gradients  near  the  droplet  interface  are  expected  to  be  very  steep,  the  grids 
are  clustered  around  it,  which  is  determined  by  a  grid  generation  function  given  later  in 
Eq.  (2.83).  The  grids  are  shown  in  Fig.  2.4. 


Figure  2.4  Illustration  of  numerical  grids  generated 
When  the  droplet  vaporizes,  its  radius  generally  decreases  with  time.  In  order  to 
implement  the  interfacial  boundary  conditions,  we  need  to  track  the  interface  during  the 
entire  computational  time.  In  this  work,  all  the  control  surfaces  are  kept  moving  with  the 
droplet  surface  during  the  vaporization  process.  The  movement  of  these  surfaces  is 


defined  by  the  following  equation: 


(2.83) 
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where  ^  is  the  transformed  grid  coordinate  quantity  whose  semi-integer  values 
correspond  to  the  positions  of  the  control  surfaces,  r  is  the  radius  of  the  control  surface, 
which  is  a  function  of  both  position  and  time.  R  is  the  radius  of  the  droplet  interface, 
which  only  varies  with  time.  /(^)  is  the  function  used  to  generate  the  grid.  Using  Eq. 
(2.83),  all  the  control  surfaces  are  moving  with  the  droplet  interface,  and  fine  grids  are 
always  kept  around  the  interface,  as  determined  by  the  function  /(^) . 

2.7.2  Independent  Variables  and  Numerical  Linearizations 

The  independent  variables  in  this  work  are  defined  as 

Q  =  ipV,pvA,peV,pYy,Vj  (2.84) 


where  V  is  the  control  volume,  A  the  area  of  the  control  surface,  and  v,  the  velocity  in 
the  r  direction.  The  volume  V  is  chosen  to  be  a  variable  since  in  the  present  moving 
boundary  problem,  it  varies  with  time.  Its  movement  is  governed  by  the  following 
equation: 


dt 


f dV  -  { w-  hdA 
%,(0 


=  0 


(2.85) 


where  w  is  the  velocity  of  control  surfaces.  For  this  one-dimensional  problem,  the 


control  volume  can  be  expressed  as 


=1®-^  t+i/2  i{-i/2= [/^f + 1  /2)  -  -inm  (2.86) 


dt  dt 


(2.87) 


where  is  the  volume  of  the  droplet. 
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Equation  (2.87)  clearly  indicates  that  all  the  control  volumes  can  be  related  to  the 
droplet  volume,  which  is  determined  by  the  interfacial  boundary  conditions. 

The  following  relation  can  be  further  developed  for  the  moving  velocity  of  a 
control  surface,  which  is  required  in  the  expression  of  the  surface  fluxes: 


w  A 


l#+l/2 


dr 

~dt 


-■ 


dt 


dV 

=  r(^  +  l/2)^ 

dt 


(2.88) 


For  application  of  the  implicit  scheme,  all  the  variables  in  the  discretized 
governing  equations  need  to  be  linearized  and  related  to  the  independent  variables  shown 
in  Eq.  (2.84).  This  linearization  processes  lead  to  the  numerical  Jacobians.  Combining 
Eqs.  (2.87)  and  (2.88),  the  flux  terms  resulting  from  grid  movement  in  the  conservation 
equations  can  be  related  to  one  independent  variable,  the  control  volume.  The 
linearization  processes  with  temperature,  pressure,  and  chemical  potential  of  each  species 
are  given  in  detail  in  the  following  subsections. 

2.7.2.1  Linearization  of  Temperature 

For  this  linearization  process,  a  thermodynamic  relation  between  temperature  and 
the  independent  variables  are  required.  First  we  begin  with  the  following  thermodynamic 
relation: 

pe  =  pe{T,p.)  (2.89) 


where  /  =  1,  •  •  • ,  . 


Its  differential  form  is 


dpe  = 


dpe 


N 


N 


dT  +  Yj^jdpj  =  pCydT  +  ^ejdpj 
JPj  7=1  7=1 


(2.90) 
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where  e.  is  the  partial  density  internal  energy  as  defined  earlier  in  this  chapter.  It  is  very 
easy  to  find 

dpeV  =Vdpe  + pedV  (2.91) 

After  substituting  Eq.  (2.90)  into  Eq.  (2.91)  and  a  little  manipulation,  the 
following  relationship  is  established: 

dpeV  =  pVC,dT  +  e.dpy  +  (pe-'^  ey.  )dV  (2.92) 

y=i  ;=i 

This  is  still  not  what  we  need  since  the  independent  variable,  pYy  in  Eq.  (2.84), 
only  ranges  from  1  to  (N-1).  Equation  (2.92)  is  further  manipulated  to  take  into  account 
of  this  difference, 

N-\  N 

dpeV  =  pVCydT  +  Cf^dpV  +  ^ (^y  "  )dpYy  +  (pe  -  ^ ^jPj )dV  (2.93) 

y=l  y=l 

Basing  on  Eq.  (2.93),  the  following  Jacobian  matrix  is  obtained: 


dQ 


_ 

pVCv 


Yj^jPj  ~  p^ 


(2.94) 


Therefore,  the  temperature  T  at  a  new  time  step  n+1  can  be  expressed  as 


2.1. 2.2  Linearization  of  Pressure 

According  to  fundamental  thermodynamic  theory,  the  following  relationship 

holds: 

P  =  p{T,Pi)  (2.96) 


(2.95) 

dQ 


Its  differential  form  can  written  as 
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dp  = 


f 


K^T'jPj 


N 


dT+'Z 


^  dp  ^ 


dpi 


(2.97) 


where  i,  j  -\,---,N .  This  relation  can  be  related  to  the  independent  variables  Q  as 


dp  = 


K^T'jPj 


N 


dT+Z 


^  dp 


1  1  ^ 

-dPiV - ZPi 

^  V  /=! 


jT,pj^i 


K^PijT,pj^i 


dV  (2.98) 


The  third  term  in  Eq.  (2.98)  can  be  further  simplified  based  on 

f  3.,  A 


i=\ 


K^Pi  jT,Pj^i 


^P\ 


y^Pj 


(2.99) 


r,K 


j 


Substituting  Eq.  (2.99)  into  Eq.  (2.98),  and  changing  the  range  of  i  in  the  second 
term  in  Eq.  (2.98)  as  1  to  (N-1),  we  have 


dp  = 


yp^ 
y^jPi 


N-\ 


dT+  Z 


^  dp  ^ 


t=l  jT,Pj^i 


^dpYiV* 


^  dp  ^ 


pPN  JT.PjtN 


dp  ^ 

ydp  jty 


dV 


(2.100) 


We  then  need  to  do  a  little  manipulation  about  the  third  term  in  Eq.  (2.100) 
according  to 


N-l  N-\ 

dpj^V  =  d{p  -  ZPiW  =  dpV  -  ZdPiV 

(=1  i=l 


(2.101) 


Substituting  Eq.  (2.101)  into  Eq.  (2.100), 
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dp  = 


N-\ 

,  dT+  V 


JT,pj^T^ 


V 


dpYiV 


+ 


dp 

dpN 


)T,pj^N 


-dpV--p\ 
V  V  ' 


y^dpjTY 


dV 


(2.102) 


Therefore, 
dp  ( dp  \  dT 


dQ 


dT 


JPj 


.  dQ 


dp 

dPN 


,0,0, 


j^N 


jT,Pj^l 


dp 

dPN 


)T,pj^f^ 


V 

(  dp  ] 

_ /I 

^  dp^ 

1 

T.PyV/V-l 

[dPN  j 

5  H 

T^Pj^N 

pPj 

(2.103) 


where  the  —  term  can  be  found  from  Eq.  (2.94). 
dQ 

2.1.23  Linearization  of  Chemical  Potential 

In  this  numerical  study,  the  chemical  potential  of  each  species  is  required  for 
treating  thermodynamic  phase  equilibrium  at  the  droplet  surface  and  calculating  the 
droplet  vaporization  rate  in  the  irreversible  thermodynamic  relationship  expressed  in  Eq. 
(2.16).  We  begin  with  the  following  thermodynamic  relation: 

Pi=Pi{T,p,Yj)  (2.104) 


where  i  =  1,  A^,  y  =  1, •  •  • , //  - 1 . 


Its  differential  form  is 


dpi  = 


d_P 

^dTj 


dT  + 


P,Yj 


V  JT,Y 


N-\Y  ^ 


dP+  X 


i=i 


dPi 

dYj 
V  •'  y 


dYi 


T.PJk^j 


(2.105) 
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where  ^  =  1,  •  •  • ,  yv  - 1 . 


We  then  need  to  express  Y.  in  terms  of  the  independent  variables  Q, 


dY.=d 


^  pYW^ 
^  J 


1  p^y 

- dpYV  --^--L-dpV 

pv  ^  >  {pvy  ^ 


1  Y 

- dpYV  --^dpV 

pV  '  pV 


(2.106) 


Substituting  Eq.  (2.106)  into  Eq.  (2.105),  we  obtain 


dPi  = 


/a,,  V 


dT  + 


P,Yj 


N- 


V  jrjj  7=1 


r  ^  \ 


dP+  X 


M 


3F,- 


dpYjV 


J  jT,P,Yk^j 


N-l  Y- 


r  ^  \ 


dYj 


dpV 


jT,P,Yk^j 


(2.107) 


From  Eq.  (2.107),  the  Jocobian  matrix  for  the  chemical  potential  is  obtained  as 


dQ 


rdpy 


j 


pJj 


dT_ 

iClQ 


dpi  ^  dp 
\  Jj^Yj 


0,0,0, 


pv 


T,p,Yki.\ 


1 

’  pV 


OF/V-l 


--S-i 


^T^P^^ki^N-l  7=1^'^ 


dYj 
V  ■'  J 


T^Pd'k^j 


(2.108) 


where  the  first  two  terms  can  be  obtained  in  Eqs.  (2.94)  and  (2.103),  respectively.  The 


terms  associated  with 


can  be  calculated 


by  using  Eq.  (2.43).  All  the 


expressions  used  in  the  linearization  processes  are  presented  in  Appendix  A  in  detail. 
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2.7.3  Solution  Procedures 

In  this  work,  a  new  method  based  on  the  irreversible  thermodynamic  phenomenon 
for  treating  interfacial  boundary  conditions  is  utilized,  as  shown  in  equation  (2.16),  where 
species  vaporization  is  related  to  the  difference  of  chemical  potential.  This  equation  is 
used  in  the  species  conservation  equation  for  the  two  interfacial  grid  cells,  and  the 
vaporization  constant  /f„ap,,  is  taken  to  be  a  very  large  number  (~10^).  Therefore,  the 

chemical  potential  difference  at  the  interface  can  be  kept  as  small  as  needed,  which  leads 
to  automatic  satisfaction  of  the  thermodynamic  phase  equilibrium  conditions  regarding 
chemical  potentials.  This  method  is  possible  because  of  the  implicit  numerical  algorithm 
utilized. 

Similarly,  when  the  heat  flux  at  the  droplet  interface  is  computed  in  the  energy 
conservation  equation,  the  thermal  conductivity  at  the  droplet  surface  will  also  be  taken 
as  a  large  number,  leading  to  the  equality  of  temperature  at  both  sides  of  the  droplet. 
Based  on  the  assumption  that  the  pressure  is  a  constant  in  the  entire  flow  field,  no  special 
effort  is  necessary  for  its  treatment.  Following  this  strategy,  numerical  iteration  is  no 
longer  required  for  solving  the  interfacial  boundary  equations. 

An  upwind  scheme  is  applied  in  the  numerical  treatment  of  the  convection  terms 
in  order  to  avoid  numerical  oscillations,  especially  at  the  droplet  surface,  where  a  large 
density  gradient  may  exist. 

Using  the  linearization  processes  described  earlier,  a  set  of  descretized  equations 
is  obtained  with  the  independent  variables  shown  in  Eq.  (2.84)  as  the  only  unknowns.  For 
the  implicit  numerical  treatment,  an  efficient  numerical  solver  is  employed  to  solve  the 
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resulting  block  tri-diagonal  matrix  (Douglas  and  Gunn  1964,  Briley  and  McDonald 
1977). 

During  the  computational  process,  when  the  relative  difference  between  species 
mass  fractions  in  the  two  interfacial  cells  becomes  sufficiently  small  (less  than  5*10'^),  a 
assumption  is  made  that  the  critical  mixing  state  is  reached.  Afterwards,  the  droplet 
becomes  a  pocket  of  dense  fluid,  and  the  flow  field  is,  therefore,  in  a  single  phase.  The 
droplet  interface  is  defined  as  an  isotherm  line  at  the  critical  mixing  temperature. 


Chapter  3 


DROPLET  CLUSTER  BEHAVIOR 

3.1  Problem  Description 

It  has  been  commonly  accepted  that  droplets  evaporate  and  bum  in  a  group  in 
practical  spray  combustion  devices.  Therefore,  there  is  increasing  interests  in  studying 
droplet  cluster  behavior  in  the  combustion  community.  This  problem  is  so  complicated 
that  some  necessary  assumptions  are  inevitable.  In  the  literature,  there  are  two  numerical 
models  for  studying  droplet  interactions  in  a  cluster.  One  is  the  “droplet-in-bubble” 
model  (Tishkoff  1979);  another  is  the  “sphere  of  influence”  concept  (Bellan  and  Cuffel 
1983).  The  two  models  are  very  similar  in  the  fact  that  they  both  take  the  following 
assumptions:  1)  The  monosize  spherical  droplets  distribute  uniformly  in  a  droplet  cluster; 
2)  There  is  a  finite  quiescent  surrounding  environment  for  each  droplet  (a  “bubble”  or 
“sphere  of  influence”),  within  which  a  numerical  method  developed  for  an  isolated 
droplet  vaporization  in  a  quiescent  environment  is  still  valid.  The  only  difference 
between  the  two  methods  is  the  treatment  of  the  interstitial  regions,  which  serve  as  the 
outside  boundary  conditions  for  those  finite  quiescent  bubbles  or  spheres.  In  the  “droplet- 
in-bubble”  model,  the  interstitial  regions  are  incorporated  into  each  bubble,  while  they 
are  treated  separately  in  the  “sphere  of  influence”  concept,  which  assumes  uniform 
variables  in  the  interstitial  regions.  These  variables  are  then  solved  using  the  global 
conservation  equations.  The  assumption  of  uniform  variables  in  the  interstitial  regions 
means  that  the  effects  of  energy  and  mass  exchanges  between  the  droplet  cluster  and  its 
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outside  environment  instantly  spread  over  the  entire  interstitial  regions.  This  assumption 
is  strictly  valid  at  the  outer  region  of  a  droplet  cluster. 

Although  the  assumptions  in  the  two  numerical  models  can  not  be  totally  justified 
for  the  practical  situations,  they  serve  to  simplify  the  droplet  cluster  problem  and  make 
numerical  solutions  possible.  Results  can  enhance  our  deep  understanding  of  the  physics 
in  droplet  interactions. 

Recently,  a  few  researchers  applied  the  “sphere  of  influence”  concept  to  study 
droplet  interactions  at  high  pressures  (Jiang  and  Chiang  1994,  1996,  Harstad  and  Bellan 
1998).  Jiang  and  Chiang  (1994,  1996)  considered  n-pentane  droplets  vaporization  in 
nitrogen,  and  Harstad  and  Bellan  (1998)  investigated  oxygen  droplet  vaporization  in 
supercritical  hydrogen  environments.  Some  interesting  results,  which  show  the 
characteristics  of  droplet  interactions  distinct  from  that  of  an  isolated  droplet,  were 
presented.  However,  considering  the  complexity  of  the  problem  and  the  simplifications 
accepted,  these  results  are  not  complete  and  conclusive.  Further  studies  based  on  different 
views  of  the  problem  are  absolutely  necessary  and  are  expected  to  complement  the 
current  results. 

In  this  chapter,  a  cluster  of  liquid  oxygen  (LOX)  droplets  vaporization  in 
hydrogen  environments  at  both  sub-  and  super-critical  conditions  is  numerically  studied. 
Figure  3.1  shows  the  schematics  of  the  problem,  where  monosized  spherical  droplets  are 
uniformly  distributed  in  a  droplet  cluster.  Since  the  assumption  of  uniform  variables  in 
the  interstitial  regions  in  the  early  works  is  valid  in  the  outer  region  of  a  droplet  cluster, 
this  study  is  focused  on  droplet  interactions  in  the  inner  region,  as  shown  in  Fig.  3.1a.  In 
this  study,  an  isolated  impermeable  bubble  is  assumed  to  surround  each  droplet  in  the 
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inner  region  of  a  droplet  cluster,  and  there  is  no  heat  and  mass  transfer  across  the  bubble 
surface.  This  is  a  reasonable  assumption  because  it  really  takes  very  long  time,  compared 
with  the  droplet  lifetime,  for  the  outside  thermal  wave  to  penetrate  into  the  droplet 
cluster.  Mass  transfer  is  also  weak  when  far  from  the  outside  boundary  (Annamalai  and 
Ryan  1992).  For  example,  with  an  oxygen  droplet  cluster  of  1cm  in  diameter  vaporization 
in  a  hydrogen  environment,  the  thermal  penetration  time  is  of  the  order  of  1000  ms,  while 
the  lifetime  of  a  droplet  at  100  pm  in  diameter  is  only  around  10  ms  at  high  pressures. 
Just  like  Tishkoff  s  model,  the  interstitial  regions  are  uniformly  incorporated  into  each 
bubble,  as  shown  in  Fig.  3.1b,  where  the  bubble  size  is  a  little  larger  than  that  of  a 
“sphere  of  influence”,  whose  radius  is  the  half  distance  between  two  neighboring  droplets 
(Bellan  and  Cuffel  1983).  Within  each  bubble,  the  environment  is  considered  as 
quiescent.  The  bubble  size  represents  the  intensity  of  droplet  cluster  effect.  For  example, 
a  small  bubble  size  means  more  droplets  in  a  droplet  cluster,  which  leads  to  stronger 
droplet  interactions.  Results  from  this  study  are  expected  to  complement  those  from  the 
early  research  works  (Jiang  and  Chiang  1996,  Harstad  and  Bellan  1998). 


Fig.  3.1  (a)  Schematics  of  a  droplet  cluster  and  the  inner  region; 
(b)  Illustration  of  the  numerical  models 
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3.2  Numerical  Treatment 

Droplet  cluster  behavior  with  liquid  oxygen  (LOX)  droplets  vaporization  in 
hydrogen  environments  at  both  sub-  and  super-critical  conditions  are  systematically 
investigated  using  the  numerical  algorithm  developed  in  Chapter  2.  The  high-pressure 
effects,  including  thermodynamic  nonideality  and  transport  property  anomaly,  are  fully 
considered.  The  computational  domain  is  shown  in  Fig.  3.2,  where  a  droplet  is 
evaporating  in  an  isolated  impermeable  bubble. 


Fig.  3.2  Single  droplet  vaporization  in  an  isolated  impermeable  bubble. 

The  bubble  size,  which  indicates  the  intensity  of  droplet  interactions,  is  initially 
prescribed,  but  it  will  change  during  the  vaporization  process,  depending  on  the 
conservation  equations  of  mass,  energy,  and  the  equation  of  state.  The  boundary 
conditions  for  this  elastic,  isolated,  and  impermeable  bubble  are  defined  as 

vr  =  o 

vy;.  =0  (3.1) 

V  =  iv 

As  described  in  Chapter  2,  all  the  control  surfaces  are  moving  with  the  droplet 
surface  during  numerical  computation,  and  the  movements  are  generally  inward  because 
of  droplet  vaporization.  Since  the  movement  of  the  control  surface  next  to  the  bubble 
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surface  is  faster  than  that  of  the  bubble  surface,  grids  become  increasingly  coarser  in  that 
region.  New  grids  have  to  be  added  in  order  to  ensure  computational  accuracy,  and  a 
simple  linear  interpolation  method,  which  proves  to  be  efficient  and  accurate,  is  adopted 
for  calculating  the  new  variables. 

3.3  Results  and  Discussions 

Droplet  cluster  behavior  at  both  sub-  and  super-critical  conditions  is  studied 
numerically  with  liquid  oxygen  (LOX)  droplet  vaporization  within  a  bubble  of  gaseous 
hydrogen,  as  shown  in  Fig.  3.2.  The  initial  temperature  and  diameter  of  the  liquid  oxygen 
(LOX)  droplet  are  lOOK  and  100|xm,  respectively,  while  the  initial  ambient  temperature 
of  hydrogen  ranges  from  700K  to  1500K.  Ambient  pressures  at  10,  30  and  100  atms, 
among  which  the  first  two  cases  represent  subcritical  states  while  the  last  one  is  at 
supercritical  condition  (see  Fig.  2.2),  are  considered.  Computation  goes  on  until  all  the 
droplet  mass  is  evaporated.  Complete  vaporization  is  considered  to  be  impossible  once 
the  saturation  conditions  are  reached.  These  are  defined  as  mass  fraction  ratio  (Y02,  inr 
Yo2.  bub)/  Yo2,  im  Icss  than  2%  or  the  temperature  at  the  bubble  surface  minus  the 
temperature  at  the  droplet  interface  lower  than  5K.  In  the  above  definitions,  Y02,  im  is  the 
mass  fraction  of  oxygen  at  the  droplet  interface  on  the  gaseous  side,  and  Y02,  bub  is  the 
mass  fraction  of  oxygen  at  the  bubble  surface.  This  research  work  is  focused  on  the 
effects  of  pressure  and  temperature  on  droplet  interactions. 

For  the  purpose  of  comparison,  droplet  lifetimes  of  an  isolated  oxygen  droplet 
vaporization  in  quiescent  hydrogen  environments  at  different  pressures  and  temperatures 
are  first  presented.  Figure  3.3  clearly  indicates  that  droplet  lifetimes  decrease 
monotonically  with  the  increase  of  pressures  or  temperatures. 
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Fig.  3.3  Droplet  lifetimes  with  an  isolated  oxygen  droplet  vaporization  in  hydrogen 
environments  at  different  pressures  and  temperatures 

In  this  study,  a  droplet  is  assumed  to  vaporize  as  an  isolated  droplet  when  the 
initial  radius  of  the  outside  bubble  reaches  100  times  of  that  of  the  droplet.  The  validity 
of  the  assumption  is  demonstrated  in  Figs.  3.4  and  3.5,  where  the  initial  temperature  and 
pressure  are  1000  K  and  30  atm,  respectively.  Figure  3.4  illustrates  the  temporal 
variations  of  temperature  outside  the  droplet  at  three  instants  when  20%,  40  %,  and  80% 
of  the  droplet  mass  are  vaporized.  Results  indicate  that  the  temperature  gradients 
approach  zero  within  a  computational  domain  whose  radius  is  less  than  50  times  of  the 
initial  droplet  radius.  The  same  conclusion  can  be  drawn  from  the  distributions  of  oxygen 
mass  fraction,  as  shown  in  Fig.  3.5 
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Fig.  3.4  Temporal  variations  of  temperature  at  an  ambient  temperature  of  1000  K  and  a 
pressure  of  30  atm 


Fig.  3.5  Temporal  variations  of  oxygen  mass  fraction  at  an  ambient  temperature  of  1000 
K  and  a  pressure  of  30  atm. 


The  effects  of  droplet  interactions  on  droplet  lifetimes  at  two  different  pressures, 
namely  30  and  100  atms,  are  presented  in  Fig.  3.6  and  3.7,  respectively,  where  the 
temporal  variations  of  nondimensional  droplet  diameter,  D^/Dq^,  are  illustrated.  In  each 
case,  droplet  vaporization  rates  at  different  ratios  of  initial  bubble  to  droplet  radius,  which 
is  called  interactive  radius  ratios,  (Rbut/Rdlo,  are  compared.  Fig.  3.6  indicates  that  at  a 
pressure  of  30atm,  the  effect  of  droplet  interactions  on  droplet  vaporization  rate  is  weak, 
once  the  interactive  radius  ratio  reaches  10.  A  droplet  behaves  exactly  like  an  isolated 
droplet  and  follows  the  traditional  law.  With  the  decrease  of  the  interactive  radius 
ratio,  droplet  vaporization  rates  change  dramatically.  At  a  lower  interactive  radius  ratio, 
or  in  a  denser  droplet  cluster,  the  droplet  vaporization  rate  follows  that  of  an  isolated 
droplet  only  at  the  short  beginning  period.  Afterwards,  the  vaporization  rate  slows  down. 
As  illustrated  in  Fig.  3.6,  there  is  a  significant  decrease  of  vaporization  rate  at  an 
interactive  radius  ratio  of  5,  as  time  increases. 

The  pressure  effect  on  droplet  cluster  behavior  is  further  illustrated  in  Fig.  3.7, 
where  pressure  is  increased  to  lOOatm.  Figure  3.7  shows  very  weak  droplet  interactions 
for  all  the  interactive  radius  ratios  (e.g.,  from  5-10). 

Further  studies  of  pressure  effects  on  droplet  lifetimes  are  thoroughly  conducted, 
and  the  results  are  presented  in  Figs.  3.8  and  3.9,  which  correspond  to  two  different 
ambient  temperatures  at  lOOOK  and  1500K,  respectively.  A  dimensionless  droplet 
vaporization  time  is  defined  as 

_  vaporation  time  of  an  isolated  droplet 
droplet  lifetime  at  a  interactive  radius  ratio 


(3.2) 
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Fig.  3.6  Transient  variations  of  dimensionless  LOX  droplet  diameter  squares  at  different 
interactive  radius  ratios  with  a  pressure  of  30atm  and  an  ambient  temperature  of  lOOOK 


Fig.  3.7  Transient  variations  of  dimensionless  LOX  droplet  diameter  squares  at  different 
interactive  radius  ratios  with  a  pressure  of  lOOatm  and  an  ambient  temperature  of  lOOOK 
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Results  in  these  two  figures  indicate  similar  trends  at  both  temperatures,  which 
strongly  suggest  that  the  ambient  pressure  presents  a  significant  effect  on  droplet 
interactions.  At  a  lower  pressure,  droplet  interactions  become  stronger.  For  example,  at  a 
given  Rbub/Ro,  is  smaller  at  a  lower  pressure,  since  the  droplet  lifetime  is  significantly 
prolonged  compared  with  that  of  an  isolated  droplet.  In  addition,  the  critical  interactive 
radius  ratio,  (Rbub/Ro)crit,  at  which  the  complete  droplet  vaporization  cannot  be  reached, 
decreases  significantly  with  increasing  pressure.  At  both  ambient  temperatures,  the 
critical  radius  ratios  at  lOOatm  are  only  around  2.5,  but  the  corresponding  values  at  lOatm 
are  larger  than  6. 

Increase  of  ambient  temperature  produces  slightly  stronger  droplet  interactions,  as 
shown  in  Figs.  3.10  and  3.11.  At  a  given  pressure,  although  the  lifetime  of  an  isolated 
droplet  is  shorter  at  a  higher  ambient  temperature,  as  shown  in  Fig.  3.3,  the  droplet  in  the 
inner  region  of  a  dense  droplet  cluster  experiences  more  lifetime  prolongation,  which 
could  compromise  the  benefit  resulting  from  high  ambient  temperature.  The  effect  of 
ambient  temperature  on  droplet  lifetime  is,  however,  significantly  weaker  than  that  of 
pressure,  as  clearly  supported  by  the  results  in  Figs.  3.8-3.11.  With  the  increase  of 
pressure,  the  effect  of  ambient  temperature  is  further  decreased.  As  illustrated  in  Fig. 
3.11,  the  differences  between  the  results  at  700K  and  lOOOK  is  negligible. 
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Fig.  3.8  Effect  of  droplet  interactions  on  LOX  droplet  vaporization  times  at  lOOOK. 


Fig.  3.9  Effect  of  droplet  interactions  on  LOX  droplet  vaporization  times  at  1500K. 
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Fig.  3.10  Effect  of  droplet  interactions  on  LOX  droplet  vaporization  times  at  a  pressure  of 
30atm. 


Fig.  3.11  Effect  of  droplet  interactions  on  LOX  droplet  vaporization  times  at  a  pressure  of 
lOOatm. 
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The  effect  of  droplet  interactions  on  the  transient  development  of  droplet  surface 
temperature  at  different  pressures  is  presented  in  Fig.  3.12.  At  a  pressure  of  lOOatm,  the 
droplet  interactions  produce  no  effect  on  the  attainability  of  the  critical  mixing  state, 
which  is  reached  after  a  very  short  time  for  both  cases  of  an  isolated  droplet  and  a  droplet 
at  a  interactive  radius  ratio  of  3.  At  a  pressure  of  30atm,  the  temporal  variations  of 
surface  temperature  are  presented  for  both  an  isolated  droplet  and  a  droplet  at  an 
interactive  radius  ratio  of  5.  Surface  temperatures  for  both  cases  increase  continuously, 
indicating  unsteady  droplet  vaporization  processes.  The  increase  rate  of  the  interacting 
droplet  is  slower,  but  eventually  the  surface  temperature  reaches  the  same  value  as  that  of 
the  isolated  droplet  at  the  end  of  the  vaporization  process.  The  situations  at  a  pressure  of 
lOatm  are  similar  to  those  at  30atm,  except  that  the  transient  variations  of  surface 
temperatures  are  weaker,  and  quasi-steady  state  surface  temperatures  are  attained. 

The  effect  of  droplet  interactions  on  temperature  distributions  outside  of  the 
droplet  is  compared  in  Fig.  3.13,  where  results  are  presented  at  three  instants,  when  20%, 
40%,  and  80%  of  droplet  mass  are  evaporated.  At  an  interactive  radius  ratio  of  5, 
temperature  outside  of  the  droplet  decreases  faster  than  that  of  an  isolated  droplet.  At  all 
three  times,  the  temperature  gradients  at  droplet  surface  are  smaller  for  the  interacting 
droplet.  This  conclusion  is  consistent  with  the  result  of  the  slower  increase  of  droplet 
surface  temperature.  Temperature  decreases  uniformly  within  the  bubble  due  to  the  large 
thermal  diffusivity  of  hydrogen. 
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Figure  3.14  presents  the  temperature  distributions  inside  the  droplet  at  the  same 
three  instants  as  those  in  Fig.  3.13.  The  transient  heating  of  the  liquid  phase  is  clearly 
illustrated  for  both  the  isolated  and  the  interacting  droplets.  With  the  increase  of  time, 
differences  of  temperature  distributions  of  the  two  cases  increase  dramatically.  For  a 
droplet  in  a  dense  droplet  cluster,  more  energy,  which  is  transferred  from  the  outside 
high-temperature  environment  into  the  droplet,  is  spent  for  heating  the  liquid  phase.  This 
leads  to  less  temperature  difference  inside  the  droplet. 

Figure  3.15  compares  the  distributions  of  oxygen  mass  fractions  outside  the 
droplet  for  the  same  two  cases  as  those  in  Fig.  3.14.  With  more  oxygen  mass  evaporated, 
the  oxygen  mass  fraction  increases  dramatically  for  the  interacting  droplet  case.  At  the 
instant  when  80%  of  droplet  mass  is  evaporated,  the  differences  of  the  oxygen  mass 
fractions  at  droplet  and  bubble  surfaces  become  extremely  close  to  the  saturation 
condition.  With  the  decrease  of  mass  fraction  gradient  at  the  droplet  surface  as  time 
increases,  the  slow  mass  transfer  rate  at  the  droplet  surface  results  in  more  oxygen  mass 
accumulation.  Because  the  thermal  conductivity  of  oxygen  is  much  smaller  than  that  of 
hydrogen,  more  oxygen  mole  fraction  at  the  droplet  surface  also  means  a  slower  heat 
transfer  rate.  Some  properties  of  pure  oxygen,  pure  hydrogen,  and  oxygen  and  hydrogen 
mixture  at  equal  molar  fractions  are  given  in  Table  3.1,  where  a  temperature  of  200K  is 
chosen  as  the  characteristic  temperature  near  the  droplet  surface. 


74 


Table  3.1a  Pressure  Effects  on  the  Properties  of  Pure  Oxygen  at  200K 

P  (atm) 

p  (kg/m^) 

Cp  (J/kg/K) 

1  (W/m/K) 

Dij  (m'/s) 

10 

20.12 

958.8 

0.201E-01 

0.393E-05 

30 

64.53 

1079.9 

0.222E-01 

0.123E-05 

100 

274.43 

1874.8 

0.331E-01 

0.284E-06 

200 

550.42 

1996.5 

0.547E-01 

0.136E-06 

Table  3.1b  Pressure  Effects  on  the  Properties  of  Pure  Hydrogen  at  200K 

P  (atm) 

p  (kg/m^) 

Cp  (J/kg/K) 

1  (W/m/K) 

Dij  (mVs) 

10 

1.22 

13901.5 

0.138 

0.408E-05 

30 

3.58 

14020.3 

0.141 

0.138E-05 

100 

11.18 

14350.8 

0.151 

0.430E-06 

200 

20.41 

14660.6 

0.164 

0.215E-06 

Table  3.1c  Pressure  Effects  on  the  Properties  of  Oxygen  and  Hydrogen 
_ at  Equal  Molar  Fractions  at  200K  _ 


P  (atm) 

p  (kg/m^) 

Cp  (J/kg/K) 

1  (W/m/K) 

Dij  (m'/s) 

10 

10.40 

17li2 

0.401E-01 

0.402E-05 

30 

31.33 

1784.2 

0.431E-01 

0.133E-05 

100 

103.8 

2006.0 

0.538E-01 

0.388E-06 

200 

195.6 

2192.2 

0.696E-01 

0.198E-06 

Further  comparisons  of  temperature  and  oxygen  mass  fraction  distributions 


outside  a  droplet  at  different  temperatures  are  presented  in  Figs.  3.16-3.17  for  an 
interacting  droplet  with  an  interactive  radius  ratio  of  5,  when  80%  of  the  droplet  mass  is 
evaporated.  Fig.  3.16  shows  different  temperature  gradients  at  the  droplet  surface,  where 
a  larger  temperature  gradient  maintains  for  an  initially  higher  ambient  temperature.  This 


is  inconsistent  with  the  former  conclusion  of  slightly  stronger  droplet  interaction  at  a 
higher  ambient  temperature.  In  contrast,  Fig.  3.17  indicates  an  opposite  trend  for  oxygen 
mass  fraction  distributions,  where  the  mass  fraction  gradient  at  the  droplet  surface  is 
smaller  for  a  higher  ambient  temperature.  This  results  in  more  oxygen  accumulation. 


Therefore,  the  slightly  stronger  droplet  interactions  at  a  higher  ambient  temperature  is 
solely  due  to  oxygen  mass  accumulation  at  the  droplet  surface,  which  decreases  the  heat 
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Fig.  3.12  Transient  variations  of  droplet  surface  temperatures  at  different  pressures  for 
both  an  isolated  and  an  interacting  LOX  droplet. 


Radius  Ratio  (r/R^) 


Fig.  3.13  Temperature  distributions  outside  the  droplets  at  a  pressure  of  30atm  and  a 
ambient  temperature  of  lOOOK  for  both  an  isolated  and  an  interacting  LOX  droplet  (1: 
20%,  2:  40%,  3:  80%  droplet  mass  vaporized) 
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Radius  Ratio  (r/R^) 


Fig.  3.14  Temperature  distributions  inside  the  droplets  at  a  pressure  of  30atm  and  a 
ambient  temperature  of  lOOOK  for  both  an  isolated  and  an  interacting  LOX  droplet  (1: 
20%,  2;  40%,  3:  80%  droplet  mass  vaporized) 


Fig.  3.15  Oxygen  mass  fraction  distributions  outside  the  droplets  at  a  pressure  of  30atm 
and  a  ambient  temperature  of  lOOOK  for  both  an  isolated  and  an  interacting  LOX  droplet 
(1:  20%,  2;  40%,  3;  80%  droplet  mass  vaporized) 
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transfer  rate  into  the  liquid  droplet  by  decreasing  the  thermal  conductivity  of  the  mixture. 
This  phenomenon  is  so  strong  that  it  totally  counteracted  the  effect  of  the  high 
temperature  gradient. 

Pressure  effects  on  temperature  and  oxygen  mass  fraction  distributions  are 
presented  in  Fig.  3.18-3.19  for  an  interacting  droplet  at  an  interactive  radius  ratio  of  5. 
Both  the  temperature  and  oxygen  mass  fraction  gradients  at  the  droplet  surfaces  are  in 
favor  of  the  high-pressure  case.  This  result  clearly  explains  the  weaker  droplet 
interactions  at  a  higher  pressure. 

Summarizing  the  early  results,  droplet  interactions  in  the  inner  region  of  a  droplet 
cluster  results  in  significantly  prolongation  of  droplet  lifetimes,  especially  at  a  small 
interactive  radius  ratio.  Lifetime  prolongation  is  mainly  caused  by  the  following  reasons: 
1)  The  rapid  decrease  of  the  bulk  temperature  in  the  bubble  results  in  decreased 
temperature  gradients  at  the  droplet  surface  during  the  droplet  vaporization  process,  and 
consequently  the  slow  heat  transfer  rate  to  the  droplet;  2)  More  energy  is  transferred 
inside  the  droplet  for  liquid  heating  for  an  interacting  droplet;  and  3)  At  the  later  stage  of 
the  droplet  vaporization  process,  low  oxygen  mass  fraction  gradient  at  the  droplet  surface 
leads  to  oxygen  accumulation  and  consequently  low  thermal  conductivity  of  the  mixture. 
The  last  reason  is  the  major  one  for  low  droplet  vaporization  rate,  or  strong  droplet 
interactions,  at  the  later  stage  of  the  vaporization  process. 

In  the  inner  region  of  a  droplet  cluster,  there  is  no  outside  energy  and  hydrogen 
supply  during  the  entire  droplet  vaporization  process.  This  conclusion  leads  to  the 
assumption  of  an  isolated  and  impermeable  bubble.  The  decrease  of  bulk  temperature  and 
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r/R, 

Fig  3.16  Effects  of  ambient  temperatures  on  temperature  distributions  outside  an 
interacting  LOX  droplet  at  a  pressure  of  30atm  when  80%  droplet  mass  vaporized. 


r/R, 


Fig  3.17  Effects  of  ambient  temperatures  on  oxygen  mass  fraction  distributions  outside 
an  interacting  LOX  droplet  at  a  pressure  of  30atm  when  80%  droplet  mass  vaporized. 
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Fig  3.18  Effects  of  pressures  on  temperature  distributions  outside  an  interacting  LOX 
droplet  at  an  ambient  temperature  of  lOOOK  when  80%  droplet  mass  vaporized. 


Fig  3.19  Effects  of  pressures  on  oxygen  mass  fraction  distributions  outside  an  interacting 
LOX  droplet  at  an  ambient  temperature  of  lOOOK  when  80%  droplet  mass  vaporized. 
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the  increase  of  oxygen  mass  fraction  are  decided  by  the  initial  total  energy  and  hydrogen 
mass  quantity  within  the  bubble.  The  initial  total  energy  depends  on  initial  hydrogen 
temperature,  as  well  as  the  initial  hydrogen  mass  quantity.  At  a  given  interactive  radius 
ratio  or  bubble  size,  temperature  and  pressure  effects  on  the  hydrogen  density  play  a 
decisive  role  for  droplet  cluster  behavior.  Pressure  effect  on  the  densities  of  initial 
gaseous  hydrogen  (e.g.,  at  lOOOK)  and  liquid  oxygen  (e.g.,  at  lOOK)  are  presented  in 
Table  3.2,  where  pressure  presents  only  slight  influence  on  the  density  of  liquid  oxygen 
but  strong  effect  on  the  density  of  hydrogen.  With  the  increase  of  pressure,  the  hydrogen 
density  increases  almost  linearly.  Therefore,  at  a  higher  pressure,  more  hydrogen  mass,  as 
well  as  more  total  energy,  strongly  slow  down  the  decrease  rate  of  bulk  temperature  and 
the  increase  rates  of  oxygen  accumulation. 

The  temperature  effect  on  the  hydrogen  density  is  shown  in  Table  3.3,  where  the 
hydrogen  density  increases  with  the  decrease  of  temperature.  However,  because  of  the 
narrow  ranges  of  temperature  variations,  its  effect  on  hydrogen  density  is  much  weaker 
than  that  of  pressure,  and  the  increase  of  hydrogen  density  is  not  capable  of  changing  the 
total  energy  sufficiently,  as  indicated  in  Fig.  3.16.  However,  the  variation  of  temperature 
can  dramatically  change  the  oxygen  mass  fraction  distribution  at  the  later  stage  of  the 
vaporization  process,  and  the  oxygen  accumulation  is  the  major  reason  for  determining 
droplet  vaporization  rate.  This  provides  the  explanation  for  the  slightly  weaker  droplet 
interactions  at  a  lower  ambient  temperature. 

TaWe_3.2^  Pressure  Effect  ori  the  Density  of  Oxygen  ^pc^Hydrogen 
Pressure  (atm)  Density  of  H?  at  lOOOK  (kg/m^)  Density  of  O2  at  lOOK  (kg/m^) 

10  '  '  0.240  .  '  1087.2 . 

30  0.717  1094.8 

100  2.355  1118.0 
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Table  3.3.  Temper^ure  Effect  oj^the JDensij^  ofj^drogen 
Temperature  (K)  Density  of  Ha  at  30atm  (kgW 

700  .  .  "1.022 

1000  0.717 

1500 _ 0.^79 _ 

Finally,  the  temporal  variations  of  the  bubble  radius  at  different  conditions,  which 
represent  the  possible  expansion  or  contraction  of  the  droplet  cluster  in  its  inner  region, 
are  studied  in  detail.  The  bubble  radius  variations  are  illustrated  in  Figs.  3.20  and  3.21, 
which  correspond  to  two  different  pressures  of  30atm  and  lOOatm,  respectively.  The 
bubble  radius  variation  is  due  to  two  factors:  the  decrease  of  hydrogen  temperature  and 
the  gasification  of  oxygen  species.  The  first  factor  tends  to  decrease  the  bubble  radius, 
while  the  second  one  increases  it. 

At  a  pressure  of  30atm,  as  shown  in  Fig.  3.20,  the  bubble  radii  decrease  at  all  the 
interactive  radius  ratios  studied  herein.  This  result  indicates  that  the  inner  region  of  a 
droplet  cluster  is  contracting  during  the  entire  vaporization  process.  The  decrease  of 
bubble  radius  is  stronger  at  a  lower  interactive  radius  ratio,  or  denser  droplet  cluster. 
When  pressure  is  increased  to  lOOatm,  the  phenomenon  of  bubble  contraction  still  exists, 
but  the  contraction  becomes  much  weaker  compared  with  that  at  30atm.  Results  indicate 
that  decrease  of  hydrogen  temperature  dominates  this  contraction  process. 
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Fig.  3.20  Transient  variations  of  bubble  radii  at  different  interactive  radius  ratios  with  a 
pressure  of  30atm  and  an  ambient  temperature  of  lOOOK. 


Fig.  3.21  Transient  variations  of  bubble  radii  at  different  interactive  radius  ratios  with  a 
pressure  of  lOOatm  and  an  ambient  temperature  of  lOOOK. 
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3.4  Summary  of  Results 

Droplet  cluster  behavior  with  liquid  oxygen  (LOX)  droplets  vaporization  in 
gaseous  hydrogen  environments  at  both  sub-  and  super-critical  conditions  is  numerically 
investigated  in  this  chapter.  A  high-pressure  droplet  vaporization  model  is  applied  with 
full  account  of  thermodynamic  nonideality  and  transport  anomaly.  Droplet  interactions  in 
the  inner  region  of  a  droplet  cluster  are  modeled  by  assuming  an  isolated  and 
impermeable  bubble  surrounding  each  droplet.  This  research  work  is  focused  on  the 
effects  of  ambient  pressure  and  temperature  on  droplet  interactions.  The  following 
conclusions  are  found  important: 

1)  Pressure  has  strong  effect  on  droplet  interactions,  while  temperature  presents  only  a 
minor  effect,  especially  at  a  high  pressure.  With  the  increase  of  pressure,  the  effect  of 
droplet  interactions  on  droplet  lifetime  decreases. 

2)  Detailed  analyses  indicate  that  the  prolongation  of  droplet  lifetime  in  a  dense  droplet 
cluster  is  mainly  caused  by  three  reasons.  First,  the  decrease  of  bulk  temperature 
leads  to  a  lower  temperature  gradient  at  droplet  surface  and,  consequently,  the  lower 
heat  transfer  rate  into  the  droplet.  Second,  more  energy  is  transferred  inside  the 
droplet  for  liquid  heating  for  an  interacting  droplet.  Finally,  at  the  later  stage  of  the 
vaporization  process,  the  lower  gradient  of  oxygen  mass  fraction  at  the  droplet 
surface  causes  more  oxygen  accumulation  and  lower  thermal  conductivity  of  the 
mixture  at  that  region.  The  third  factor  dominates  at  the  later  stage  of  the  droplet 
vaporization  process. 

3)  The  effects  of  the  initial  ambient  pressure  and  temperature  on  hydrogen  density  play  a 
decisive  role  in  determining  the  initial  hydrogen  mass  quantity  and  total  energy  in  the 
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isolated  and  impermeable  bubble,  and  ultimately  decided  the  decrease  rates  of 
temperature  and  mass  fraction  gradients  at  the  droplet  surface. 

4)  The  inner  region  of  a  droplet  cluster  presents  continuous  contraction  during  the  entire 
droplet  vaporization  process,  and  the  phenomenon  is  stronger  at  a  lower  pressure. 

The  results  of  this  study  complement  those  from  the  early  researches  (Jiang  and 
Chiang  1996,  Harstad  and  Bellan  1998),  and  give  us  a  complete  physical  understanding 
of  droplet  cluster  behavior. 


Chapter  4 


UNIFIED  TREATMENT  OF  REAL  FLUID  MIXTURES  AND  THE 
APPLICATION  TO  PRECONDITIONING  SCHEME 

In  practical  spray  combustion  devices,  droplets  generally  evaporate  in  convective 
environments.  A  unified  numerical  algorithm  capable  of  treating  the  real  fluid  behavior 
ranging  from  dense  liquid  to  dilute  gas  has  to  be  developed.  All  the  high-pressure 
characteristics,  including  thermodynamic  nonideality  and  transport  anomaly,  must  be 
properly  taken  into  account.  In  addition,  the  numerical  algorithm  must  be  able  to  handle 
small-Mach  number  fluid  flows  in  the  dense  liquid  region. 

In  this  chapter,  a  unified  treatment  of  real  fluid  behavior  is  developed,  and  is 
applied  to  the  preconditioning  scheme.  The  unique  features  in  this  numerical  algorithm 
are  as  following.  First,  definitions  of  the  partial  mass  and  partial  density  properties  are 
introduced  to  take  into  account  of  the  thermodynamic  properties  in  real  fluid  mixtures. 
Second,  thermodynamic  properties  are  derived  using  fundamental  thermodynamic 
theories,  and  they  can  be  calculated  based  on  any  equation  of  state  to  naturally  take  into 
account  of  the  pressure  effect.  Finally,  the  development  is  applied  to  the  preconditioning 
scheme,  which  renders  the  numerical  algorithm  capable  of  solving  small-Mach  number 
fluid  flows.  With  the  combinations  of  the  consistent  thermodynamic  treatments  of  real 
fluid  behavior  and  the  preconditioning  scheme,  an  efficient  and  robust  numerical 
algorithm  is  developed  in  this  chapter  for  solving  all-Mach  number  fluid  flows  at  any 


fluid  state. 


4.1  Preconditioning  Scheme 

4.1.1  Conservation  Equations 
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The  transient  conservation  equations  of  mass,  momentum,  energy,  and  species  are 
compactly  expressed  as 

drQ  dr(E-B,)  MF-F,) 

3l  dr  ^  ’ 

where  the  equations  are  established  in  an  axisymmetric  coordinate  system  in  order  to  be 
consistent  with  the  computations  in  the  next  two  chapters. 

The  conserved  variable  vector  Q  is  defined  as 

Q  =  [p  pu^  pUr  pe,  pYiJ  (4.2) 

The  convective  flux  vectors  in  the  axial  and  radial  directions,  E  and  F,  take  the 
following  forms: 


- 1 

1 

1 

1 _ 

pHx  +P 

pu^Uj. 

pu^Uf. 

(4.3)  F  = 

2 

puj.  +  p 

(pet  +  p)u^ 

{pet  +  p)Ur 

pUxYi 

pUfYi 

(4.4) 


The  source  vector  S,  which  includes  the  terms  associated  with  the  axisymmetric 
geometry,  viscous  diffusion,  and  chemical  reaction,  is  expressed  in  the  following  form: 
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5  = 


dx 


(2 


-mr 

J 

4  uu^  2  du^  2  du 
P  - - -M, -r- 
3  r  3  dx  3  dr 
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(2  ^ 

-  ^ 

2] 

dx 

dr' 

J 

rcO: 


(4.5) 
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where  /  =  1,  •  •  •  ,  and  the  variable  N  is  the  total  number  of  species  involved  in  the 

problem. 

The  corresponding  viscous  flux  vectors,  E^,  and  F„,  can  be  written  as  second- 
order  terms.  In  order  to  do  this,  a  viscous  vector  Q,,  is  first  defined  as 


11 

(4.6) 

The  viscous  flux  vectors  are  then  expressed  as 

P  -  p  ^2.  p  ^2., 

OX  or 

(4.7) 

P  -  p  ^2.  p  dQ^. 

(4.8) 

where  the  viscous  coefficient  matrixes  are 
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(4.9a) 


(4.9b) 
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(4.9c) 


'0  0  0  0  0  •••  0 

0^  0  0  0  •••  0 

0  0  0  0  •••  0 
3 

^rr-  0  ^IJU^  A  •••  P^N-\)^{N-\)m~hiMD^f^ 

0  0  0  0  pDi„,  •••  0 

0  0  0  0  0  •  •  •  f^{N-\)m 

(4.9d) 

The  standard  variables  are  used  in  Eqs.  (4.6-4.9),  including  density  p ,  velocity 
components  in  x-  and  r-  directions  and  Uy ,  pressure  p ,  temperature  T  ,  the  total 
energy  per  unit  mass  ,  mass  fraction  of  species  i  Yi ,  viscosity  p  ,  conductivity  A , 
mass  diffusivity  in  the  mixture  ,  and  the  chemical  reaction  rate  of  species  i  cbi . 
There  is  only  one  special  variable,  the  partial  mass  enthalpy  of  species  i  /i,- ,  which  is 


defined  in  Chapter  2. 
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4.1.2  Preconditioning  Method 

It  is  well  known  that  time-marching  algorithms  are  extremely  efficient  for  solving 
the  transonic,  supersonic,  and  hypersonic  compressible  flow  problems,  but  they 
encounter  convergence  difficulties  for  small-Mach  number  flows.  Two  major  reasons 
account  for  this:  the  machine  round-off  errors  arising  from  the  calculation  of  the  pressure 
gradients  and  the  stiffness  of  the  eigenvalues  in  the  numerical  system. 

The  first  problem  can  be  circumvented  easily  by  introducing  pressure 
decomposition  (Choi  and  Merkle  1993,  Shuen  et  al.  1993), 

P  =  PQ  +  Pg  (4-10) 

where  pg  a  constant,  which  is  generally  chosen  as  the  major  part  of  the  pressure.  The 
gauge  pressure  pg  is  the  dynamic  part,  responsible  for  the  velocity-pressure  coupling  in 

the  momentum  equations,  and  it  is  the  driving  force  for  fluid  flows. 

After  pressure  decomposition,  the  total  pressure  in  the  momentum  equations  is 
replaced  by  the  gauge  pressure  pg  .  We  then  have  the  following  new  inviscid  flux  and 

source  terms: 
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fXUyUy^ 
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(4.13) 


4  uur  2  dur  2  du 

^  3  r  3  dx  3  dr 

-A 

dx  ^ 

rCO; 


The  new  viscous  variable  vector  is  also  redefined  as 

Q^=\pg,u„u„T,Yi]'  (4.14) 

All  the  other  terms  in  the  conservation  equations  are  not  affected  by  pressure 
decomposition. 

In  order  to  circumvent  the  difficulty  of  eigenvalue  stiffness,  we  need  to  add  in  a 

preconditioning  term  T— — ,  where  the  variable  r  represents  the  pseudo-time.  The 

d-i 

preconditioning  matrix  T,  and  the  variable  vector  Q  are  chosen  as  (Shuen  et  al.  1993) 
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(4.15) 


and, 


Q  =  \pa,u^,Ur,h,Yi^ 


(4.16) 


where  the  variables  are  enthalpy  h ,  the  total  enthalpy  per  unit  mass  /i, ,  and  a  scaling 


parameter  fd . 


The  conservation  equations  with  the  pseudo-time  derivative  terms  are  expressed 


„  drQ  drQ  drE  drF 

1  - 1 — - 1 - 1 - 

dt  dt  dr 


dx 
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dr 
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^Qv 
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=  5 


(4.17) 


Proof  of  the  appropriateness  of  the  chosen  preconditioning  term  for  a  real  fluid 
mixture  is  presented  in  the  following  sections,  where  the  scaling  factor  is  also  defined. 
4.2.  Numerical  Derivations  for  Real  Fluid  Mixtures 

All  the  relationships  required  in  the  numerical  algorithm  are  derived  using 
fundamental  thermodynamic  theories.  They  can  be  calculated  based  on  any  equation  of 
state.  The  partial  mass  and  partial  density  properties  of  a  real  fluid  mixture,  as  defined  in 
Chapter  2,  are  also  utilized  in  the  derivations. 

4.2.1  Important  Thermodynamic  Relationships 

Some  important  thermodynamic  relationships  are  needed  for  the  derivations  of 
numerical  Jacobians  and  eigenvalues  in  the  following  sections.  We  will  derive  these 
relationships  here. 

First,  a  thermodynamic  relationship  correlating  enthalpy  as  a  function  of  pressure, 
density,  and  mass  fractions  are  derived.  According  to  thermodynamics,  each  intensive 
property  will  depend  on  A"  -t- 1  other  intensive  variables  in  the  mixture.  We  begin  with 
the  following  relation: 

pe  =  pe(T,pi)  (4.18) 


where  /  =  I,---,  A  ,  and  e  is  the  internal  energy  per  unit  mass. 
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We  are  only  interested  in  its  differential  form,  which  can  be  expressed  as 


^de'^ 


dpe  =  p\ 


?T\ 

K^^Jpi  i=l 


dpe 


(4.19) 


Based  on  the  definition  of  partial  density  properties,  it  is  recognized 


= 


^  dpe 

y^^‘JT,pj^i 


which  is  the  partial  density  internal  energy  of  species  i  in  the  mixture.  The  first  derivative 
in  Eq.  (4.19)  is  the  constant  volume  heat  capacity  C^,.  Substituting  these  variables  into 


Eq.  (4.19)  leads  to  the  following  expression: 


N 


dpe  =  pCydT  +  Yj  ?idPi 
/=1 


(4.20) 


We  then  want  to  correlate  temperature  T  as  a  function  of  pressure  p  and  the 
partial  density  Pi .  According  to  thermodynamics,  the  following  relation  exists: 

P  =  P^^Pi)  (4.21) 

where  i  -  I,--- ,N  . 


Its  differential  form  is  expressed  as 
^dp^ 


dp  = 


dT+Y, 


Pi 


dpi 


(4.22) 


JT,pj^ 


In  this  equation,  the  summation  in  the  second  term  is  first  changed  to  be 
1,  •  •  •  A  - 1 .  Some  simple  manipulations  will  then  lead  to  the  following  relation: 
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(4.23) 

Substituting  Eq.  (4.23)  into  Eq.  (4.20),  the  following  expression  is  easily  derived: 
/  N  f 
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(4.24) 


The  left  hand-side  term  is  further  rewritten  as 


dpe  =  pde  -h  edp 


(4.25) 


Substituting  Eq.  (4.25)  into  Eq.  (4.24),  the  following  relation  is  established  after 
the  reorganization: 
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(4.26) 


Based  on  fundamental  thermodynamic  theories,  the  following  relation  can  be 


easily  obtained: 
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dh  =  de  +—dp --^dp 
P  p^ 


(4.27) 


Substituting  Eqs.  (4.26)  into  Eq.  (4.27),  the  following  expression  can  be  obtained 
after  some  straightforward  manipulations: 
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(4.28a) 

This  relation  will  be  used  frequently  for  the  derivation  of  the  numerical  Jacobians. 
The  following  compact  notations  are,  therefore,  defined: 


V  m 
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Equation  (4.28a)  is  finally  expressed  in  the  following  form: 


A^-1 


dh  =  Apdp  +  Apdp  +  Ay.  dY, 
/=!  ' 


(4.28d) 


(4.28e) 


Equation  (4.28e)  is  one  of  the  thermodynamic  relationships  we  are  looking  for. 


95 


Next,  another  relationship,  which  correlates  enthalpy  as  a  function  of  pressure, 
temperature,  and  mass  fractions,  will  be  developed.  We  start  with  the  following  relation: 
dpi  =  pdYi  +  Yidp  (4.29) 


Using  Eq.  (4.29),  Eq.  (4.20)  can  be  rewritten  as 
N  N 

dpe  =  pCydT  +  '^peidYj  + 

i=l  i=l 

Using  Eq.  (4.25),  Eq.  (4.30)  can  be  further  derived  as 


(4.30) 
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(4.31) 


Based  on  Eqs.  (4.25)  and  (4.29),  a  differential  relationship  correlating  density  as  a 
function  of  temperature,  pressure,  and  the  mass  fractions  can  be  expressed  as 


dp 


dp  = 


J  Pi 


N-l 

dT-  S  Pi 
/=! 


^  A 


dpi 


jT,Pj^i 


\dY; 
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(4.32) 


Substituting  Eqs.  (4.27)  and  (4.32)  into  Eq.  (4.31),  another  important  differential 
thermodynamic  relationship  is  established 


N-l 


dh  —  BfdT  +  B  pdp  +  By^  dY; 

i=\  ‘ 

where  the  coefficients  Bj,Bp,  and  By-  are  defined  as 


(4.33a) 
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By.  =  {Si  -Sf^)- 


(N 
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(4.33d) 

According  to  the  definition  in  thermodynamics,  we  recognize  that  the  coefficient 
Bj  is  equal  to  the  constant  pressure  heat  capacity  Cp  of  a  mixture, 


Cp=%=C„-- 
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(4.34) 


Finally  a  relationship  regarding  the  speed  of  sound  in  the  mixture  is  derived. 
According  to  the  definition  of  the  speed  of  sound. 


\^P^Js,Yi 


(4.35) 


Based  on  Eqs.  (4.25)  and  (4.29),  the  following  relationship  is  derived: 
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(4.36) 

Based  on  Eq.  (4.36),  the  following  expression  is  obtained  in  a  straightforward 


manner: 
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In  thermodynamics,  the  following  relationship  exists: 
^  =  s{T,p,Yi) 


(4.38) 


where  /  =  I,--- -1. 
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After  utilizing  some  fundamental  thermodynamic  relationships,  the  following 
differential  form  of  Eq.  (4.38)  can  be  obtained: 

ds 


ds  =  ^dT--\r(^]  dp+^ 
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i=l 
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jT,p,Yj^i 
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(4.39) 


Based  on  Eq.  (4.39),  the  following  expression  is  further  derived: 


s,Yi  P' 
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(4.40) 


Substituting  Eq.  (4.40)  into  Eq.  (4.37),  a  general  expression  of  the  speed  of  sound 
in  a  mixture  is  established  as 
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(4.41) 


\  jTJi 

Equations  (4.28),  (4.33),  (4.34),  and  (4.41)  are  the  important  thermodynamic 
relationships  required  in  the  numerical  derivations  below. 

4.2.2  Implicit  Numerical  Scheme  and  Jacobians 

The  preconditioned  transient  conservation  equations  are  solved  using  the  implicit 
scheme  with  the  dual  time-stepping  integration  technique.  The  discretized  equations  are 
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(4.42) 
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/S  A. 

where  T,  D,A,5,7?^^,and  Rj.,.  are  the  numerical  Jacobians  defined  in  the  linearization 
process; 


j.  dQ 
dQ 

dE 

(4.43)  A  =  ^ 

dQ 

dF 

(4.44)  B  =  ^ 

dQ 

(4.45) 

D=^ 

dQ 

(4.46) 

oQ 

(4.47)  = 

dQ 

(4.48) 

The  coefficients  aj ,  02 ,  and  03  are  chosen  based  on  the  requirement  of  the  time 


accuracy.  In  the  case  that  time  steps  are  uniform  and  the  second-order  time  accuracy  is 
required,  they  can  be  chosen  as 


fit] 


>«2 


-2,^3=  i 


(4.49) 


The  left  hand-side  of  Eq.  (4.42)  is  evaluated  in  the  pseudo-time  step  p .  The 
superscripts  n  +  l,n,and  n-lin  the  right  hand-side  represent  the  variables  in  the  new, 
current,  and  old  physical  time  steps,  respectively.  Once  steady  state  solution  is  reached 
through  pseudo-time  iterations,  the  new  solution  in  the  physical  time  step  n  +  1  is 
obtained. 

The  thermodynamic  relationships  derived  in  the  section  above  will  be  utilized  to 
find  the  numerical  Jacobians.  Only  the  Jacobian  A  is  presented  here  as  an  illustrative 
example,  all  the  other  Jacobians  are  given  in  Appendix  C. 

During  the  derivation  of  the  numerical  Jacobian  A,  a  relationship  correlating 
enthalpy  as  a  function  of  density,  pressure,  and  mass  fractions  is  obviously  required. 
Based  on  the  thermodynamic  relationships  shown  in  Eqs.  (4.28b)-(4.28e),  the  numerical 
Jacobian  A  takes  the  following  form: 
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4.2.3  Eigenvalues 

For  the  preconditioned  conservation  equations,  the  eigenvalues  in  the  x-direction 


are  obtained  in  the  matrix  F  ^A.  The  calculation  is  straightforward,  and  the  eigenvalues 


are 


Ux 

iul(l-Ay)^+AI3 

. 

1  J 

V 

ac 

,  ,  •  •  • 


(4.51) 


where 


a2=AV 

“c  1 

--A 
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(4.52) 


We  will  prove  that  the  variable  is  indeed  equal  to  the  speed  of  sound  a  in  a 


mixture.  The  variable  Ap  is  first  reorganized  as 
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JpX^^  JTJi  P  i=\  P  P 


Jpi  I  ^P  )T,Yi  P\^^  JPi  I  ^P  jTJi  [i=l  P  j 


It  is  recognized  that  the  term  in  the  brackets  is  equal  to  the  constant  heat  capacity 
C„  (e.g.,  referring  to  Eq.  (4.34)).  Therefore,  the  following  relation  is  derived; 


A  c 


(4.53b) 


Utilizing  Eq.  (4.28b),  it  then  gives 


l-A  --Cl^ 

p  isp 


(4.54) 


Substituting  Eqs.  (4.53b)  and  (4.54)  into  Eq.  (4.52),  a  final  relation  is  derived  as 


dT]  (dp 
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(4.55) 


The  eigenvalues  in  x-direction  can  now  be  expressed  as 


1  (  P)  2  P  2 
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^  V  /7^ 


(4.56) 


This  result  is  extremely  important  since  it  indicates  the  eigenvalues  of  the 
preconditioned  conservation  equations  of  real  fluid  mixtures  are  in  the  same  form  as 
those  of  the  ideal  gas  mixtures  (Shuen  et  al.  1993,  Buelow,  1995).  This  result  also 

demonstrates  that  our  choice  of  the  preconditioning  matrix  F  and  main  variables  Q  are 


appropriate.  The  derivations  are  independent  of  any  equation  of  state. 
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Based  on  the  eigenvalues  in  Eq.  (4.56),  another  important  conclusion  is  that  the 
scaling  factor  /5  can  be  chosen  exactly  as  that  derived  in  an  ideal  gas  mixture.  Since  the 
eigenvalues  remain  the  same  forms,  the  same  scaling  factor  will  definitely  work  for  both 
cases  to  overcome  the  eigenvalue  stiffness  problem. 

4.2.4  Selection  of  the  Scaling  Factor 

For  inviscid  flows,  in  order  to  get  the  well-conditioned  eigenvalues,  the  scaling 
factor  can  be  taken  as  (Choi  and  Merkle  1993,  Shuen  et  al.  1993,  Hsieh  and  Yang  1997) 


Pinv  ~ 


u\ 


2  2 
if  u  <uj^ 

if  ug  <u  <  a 

if 


(4.57) 


where  the  variable  uj^  is  a  reference  velocity.  Generally  a  very  small  number  is  chosen  to 
maintain  the  reference  Mach  number  at  the  order  of  10'^. 

Referring  to  the  eigenvalue  expressions  in  Eq.  (4.56),  we  found  that  this  choice  of 
scaling  factor  indeed  maintains  all  the  eigenvalues  at  the  same  order  for  inviscid  all-Mach 
number  flows.  Moreover,  all  the  eigenvalues  remain  as  real  numbers,  and  the 
preconditioned  conservation  equations  are  still  of  hyperbolic  nature. 

In  some  viscous  flow  regions,  since  the  flow  Reynolds  number  Re  can  become 
far  smaller  than  1,  the  viscous  diffusions  becomes  dominant  compared  with  the 
convection  terms.  Under  these  circumstances,  the  scaling  factor  must  also  be  able  to 
control  the  viscous  damping  rates.  Therefore,  an  appropriate  scaling  factor  for  Navier- 
Stokes  equations  should  be  able  to  simultaneously  control  both  the  Courant-Friedrichs- 
Lewy  (CFL)  and  the  von  Neumann  (VNN )  numbers.  This  can  be  accomplished  by 
choosing  the  scaling  factor  as  (Choi  and  Merkle  1993,  Buelow  1995) 
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yS  =  min(l,  max(y6,„^; ,  ))  (4.58) 

where  the  inviscid  scaling  factor  is  given  in  Eq.  (4.57).  The  viscous  scaling  factor 

will  have  to  consider  the  grid  aspect  ratio  effect,  which  is  an  important  factor  for 

viscous  flow  computation.  For  the  2-D  and  axisymmetric  flows,  the  viscous  scaling  factor 
is  chosen  as  (Buelow  1995) 


Pvis  ) 


(4.59a) 


where  the  two  scaling  factors  in  both  x-  and  r-  directions  are  defined  as 
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The  variables  o^,o^  are 
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The  other  variables  are  Courant-Friedrichs-Lewy  {CFL)  number,  the  von 
Neumann  {VNN )  number,  the  grid  aspect  ratio  AR ,  and  the  cell  Reynolds  numbers  in  x- 
and  r-  directions  Re^,  Re^^^,  respectively.  The  definitions  of  these  variables  are 
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VNN  =  max( 
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The  scaling  factor  of  the  preconditioned  conservation  equations,  which  is  valid  for 
both  inviscid  and  viscous  flows,  is  now  completely  established. 

4.3  Thermophysical  Properties 

The  Soave-Redlich-Kwong  (SRK)  equation  of  state,  as  presented  in  Chapter  2,  is 
utilized  for  calculating  all  the  thermodynamic  properties  with  full  account  of 
thermodynamic  nonideality.  Some  relationships  have  already  been  derived  in  Chapter  2, 
and  more  details  are  presented  in  Appendix  A. 

The  transport  properties  to  be  evaluated  herein  are  thermal  conductivity,  viscosity 
and  mass  diffusivity.  These  properties  are  evaluated  for  a  real  fluid  mixture  with  proper 
consideration  of  transport  anomaly.  To  evaluate  thermal  conductivity  and  viscosity,  the 
extended  corresponding  states  one-fluid  model  is  used  (Ely  and  Hanly  1981,  1983).  The 
binary  mass  diffusivities  at  high  pressures  are  evaluated  by  the  method  proposed  by 
Takahashi  (1974),  which  is  a  very  simple  corresponding  state  method.  The  details  can  be 
found  in  Chapter  2. 
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4.4  Numerical  Solution 

The  preconditioned  conservation  equations  are  discretized  by  central-differencing 
in  space.  With  an  oxygen  droplet  vaporization  in  supercritical  hydrogen  environments,  a 
contact  discontinuity  with  a  very  large  density  gradient  exists  between  the  oxygen  and 
hydrogen  interface.  Therefore,  artificial  viscous  dissipation  terms  are  required.  In  this 
numerical  treatment,  the  second-order  and  fourth-order  explicit  artificial  viscous  terms 
are  used.  An  implicit  second-order  artificial  viscous  term  is  also  added  in  the  left  hand- 
side  in  order  to  make  the  scheme  stable.  With  these  artificial  viscous  terms,  the 
discretized  governing  equations  become 
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(4.60) 

where  represents  the  implicit  second-order  artificial  viscous  terms,  and 


represent  the  explicit  second-order  and  fourth-order  artificial  viscous  terms,  respectively. 
The  implicit  second-order  artificial  viscous  terms  can  be  further  expressed  as 


(4.61) 


where  the  term  in  the  x-direction  is 


(4.62) 
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In  Eq.  (4.62),  A;,.,  and  V^.  are  the  standard  forward  and  backward  difference 
operators  in  x-direction,  and  is  the  viscous  coefficient,  which  will  be  defined  later  in 
this  section. 

The  term  in  the  r-direction  can  be  expressed  similarly. 

The  explicit  second-order  artificial  viscous  terms  can  be  expressed  as 


where  the  term  in  the  x-direction  is  given  as 

The  coefficient  £ix\+\/2  j  defined  as 


(4.63) 
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(4.65) 


where  =  0.5 . 

To  make  the  scheme  TVD,  the  following  flux  limiter  function,  ,  is  defined  as 


(Swanson  and  Turkel  19992,  Jorgenson  and  Turkel  1993): 
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Pi+\J  PiJ 

PiJ  Pi-\J 

)  +  (^(Pi  +  lJ  +  '^PiJ  +  A-1.;) 

(4.66) 

where  &)  =  0.01  is  chosen  in  this  research  work.  Density  is  chosen  as  the  controlling 
parameter  in  Eq.  (4.66)  since  contact  discontinuity  exists  in  this  problem. 

The  terms  relating  to  r-direction  can  be  expressed  similarly. 
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The  explicit  fourth-order  artificial  viscous  terms  can  be  expressed  as 

(4.67) 

where  the  term  in  the  x-direction  is  defined  as 

oS*  =  (4.68) 

The  coefficient  £eV i+\i2  j  defined  as 


^ex  i+\l2,j  ni^x(o, ^  ^ex  i+\l2,j) 


(4.69) 


The  parameter  is  a  constant,  which  is  taken  as  in  this  numerical 

algorithm. 

In  order  to  make  the  scheme  stable,  the  coefficient  in  the  implicit  second- 
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order  artificial  viscous  terms  is  chosen  as 
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The  central-difference  TVD  scheme  for  the  preconditioning  numerical  scheme  is 


completely  defined.  The  descretized  equations  are  solved  by  the  modified  strongly 
implicit  procedure  (MSIP)  originally  developed  by  Schneider  and  Zedan  (1981). 


Chapter  5 


SUPERCRITICAL  OXYGEN  DROPLET  VAPORIZING  IN 
CONVECTIVE  HYDROGEN  ENVIRONMENTS 

5.1  Problem  Description 

In  practical  spray  combustion  devices,  such  as  diesel  engines,  high-performance 
gas  turbine  engines,  and  liquid-propellant  rocket  engines,  droplets  are  vaporizing  in 
supercritical  convective  environments.  During  the  vaporization  process,  a  droplet  may 
deform  and  even  break  up.  Therefore,  a  practical  requirement  is  to  study  the  convection 
effect  on  supercritical  droplet  vaporization  processes. 

In  this  chapter,  an  oxygen  droplet  vaporizing  in  convective  hydrogen 
environments  at  supercritical  conditions  is  numerically  studied.  The  oxygen  droplet  is 
assumed  to  reach  its  critical  state  immediately  after  being  injected  into  the  combustor. 
Therefore,  it  becomes  a  problem  of  supercritical  fluid  flow  with  heat  and  mass  transfer. 
To  numerically  analyze  this  problem,  several  difficulties  have  to  be  circumvented.  First, 
the  oxygen  droplet  is  injected  into  the  combustion  chamber  at  an  approximate 
temperature  of  lOOK,  which  is  lower  than  the  critical  temperature  of  oxygen  (e.g., 
154.6K).  Therefore,  at  supercritical  conditions,  the  oxygen  droplet  and  the  oxygen- 
hydrogen  mixture  in  the  region  near  the  droplet  can  not  be  assumed  as  an  ideal  gas 
mixture.  A  consistent  and  efficient  numerical  algorithm  has  to  be  developed  to  handle  the 
real  fluid  behavior.  Second,  in  a  droplet  vaporization  and  combustion  problem,  density 
variation  is  significant  due  to  the  change  of  mass  fractions  of  the  species.  The  numerical 
algorithm  for  solving  compressible  fluid  flows  is  a  natural  choice.  However,  in  the  dense 
oxygen  droplet  region,  the  flow  velocity  is  generally  very  small.  A  numerical  method 
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capable  of  solving  small-Mach  number  flows  is  a  prerequisite  for  this  study.  Third,  for 
this  supercritical  droplet  vaporization  case,  a  very  large  density  gradient  exists  between 
the  dense  oxygen  and  the  surrounding  hydrogen.  For  example,  the  oxygen  density  is 
about  1200kg/m^  at  lOOatm  and  lOOK,  while  the  hydrogen  density  is  approximately 
2.4kg/m^  at  the  same  pressure  and  lOOOK.  In  this  case,  the  density  ratio  reaches  500  to  1. 
An  appropriate  numerical  method  should  be  utilized  to  treat  the  large  density  gradient. 
Finally,  accurate  evaluations  of  thermophysical  properties  always  present  a  great 
challenge  for  solving  supercritical  droplet  vaporization  problems. 

A  supercritical  droplet  vaporization  in  a  convective  environment  is  illustrated  in 
Fig.  5.1.  The  following  assumptions  are  taken  in  this  study: 

1 .  The  flow  is  laminar  and  axisymmetric; 

2.  The  oxygen  droplet  reaches  the  critical  point  immediately  after  being  injected  into  the 
combustion  chamber.  This  is  an  approximation  for  most  of  the  hydrocarbon  fuels,  but 
it  is  true  for  this  oxygen-hydrogen  case  (Yang  et  al.  1994,  Lafon,  1995,  Handenwang 
et  al.  1996). 


Supercritical  Stream 


Fig.  5.1  Schematic  of  a  supercritical  droplet  in  a  convective  stream. 
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5.2  Numerical  Treatment 

A  preconditioning  numerical  algorithm  with  a  unified  treatment  of  real  fluid 
mixture,  which  is  developed  in  Chapter  4,  is  utilized  to  solve  the  current  supercritical 
droplet  vaporization  problem.  A  central-difference  TVD  method  is  applied  to  handle  the 
large  density  gradient  (Swanson  and  Turkel  1992,  Jorgenson  and  Turkel  1993).  The  state- 
of-the-art  techniques  are  taken  for  thermophysical  property  evaluations,  as  presented  in 
Chapter  2. 

During  the  vaporization  process,  a  droplet  will  move  and  experience  strong 
deformation.  There  is  no  advantage  to  generate  grids  that  are  initially  conforming  to  the 
droplet  surface.  Therefore,  a  grid  system  shown  in  Fig.  5.2  is  used  in  this  numerical 
study.  The  mesh  in  the  region  around  the  droplet  is  refined  for  computational  accuracy. 
The  distances  of  the  outside  boundaries  from  the  droplet  are  at  least  15  times  of  the 
droplet  radius  to  avoid  the  far  field  effect. 


X  (m) 


Fig.  5.2  Part  of  the  computational  grids 
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5.3  Results  and  Discussions 

An  oxygen  droplet  vaporizing  in  convective  hydrogen  environments  at 
supercritical  conditions  is  shown  in  Fig.  5.1.  In  this  research  work,  the  droplet  is  50pm  in 
diameter.  Initially,  the  droplet  is  at  a  uniform  temperature  of  lOOK,  while  the  inflow 
hydrogen  is  at  lOOOK.  Pressures  range  from  100  to  400  atm.  The  flows  studied  herein  are 
restricted  to  those  with  the  Reynolds  number  (Re)  less  than  240.  The  Reynolds  number 
(Re)  is  defined  based  on  the  initial  droplet  diameter  d^,  the  inflow  hydrogen  velocity  uq  , 

and  the  hydrogen  properties  )•  The  Reynolds  number  limit  is 

based  on  the  result  of  an  incompressible  viscous  fluid  flow  past  a  stationary  sphere,  in 
which  a  steady  axisymmetric  flow  can  only  be  maintained  with  the  Reynolds  number  up 
to  200  (Nakamura  1976,  Johnson  and  Patel  1999).  For  this  supercritical  droplet 
vaporization  case,  at  the  present  time,  no  experimental  or  3-D  numerical  results  are 
available  to  validate  the  upper  limit  of  the  Reynolds  number  for  the  axisymmetric 
assumption. 

For  this  supercritical  droplet  vaporization  case,  although  the  droplet  has  already 
been  assumed  to  reach  critical  mixing  state,  its  core  region  with  low  temperatures  still 
possesses  liquid-like  density  and  characteristics.  Since  the  droplet  becomes  a  dense  fluid, 
it  leaves  us  the  freedom  to  define  the  droplet  surface.  In  the  early  works  regarding  an 
isolated  droplet  vaporizing  in  a  quiescent  environment  (Yang  et  al.  1994,  Lafon,  1995), 
an  isothermal  line  at  the  critical  mixing  temperature  of  the  oxygen-hydrogen  mixture,  as 
indicated  by  the  phase  diagrams  in  Fig.  2.2,  is  generally  defined  as  the  surface  of  the 
oxygen  droplet.  Fig.  2.2  clearly  illustrates  that,  at  a  given  pressure,  a  critical  mixing  state 
is  not  only  corresponding  to  a  critical  mixing  temperature,  but  also  to  the  critical  mixing 
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mass  fractions.  Therefore,  instead  of  the  critical  mixing  temperature,  a  critical  mixing 
mass  fraction  can  also  be  used  to  define  the  droplet  surface  (Hadenwang  et  al.  1996). 
Furthermore,  a  critical  mixing  state  in  Fig.  2.2  exists  only  at  the  states  with 
thermodynamic  phase  equilibrium.  Without  phase  equilibrium,  the  critical  temperature 
and  critical  mass  fractions  are  not  located  at  the  same  position,  rendering  the  critical 
mixing  state  meaningless.  This  is  actually  the  case  after  the  droplet  reaches  the  critical 
mixing  state  as  a  result  of  different  heat  and  mass  transfer  rates.  Therefore,  the  droplet 
surface  is  defined  as  an  isothermal  line  at  the  critical  temperature  of  oxygen  (154.6K)  in 
this  work,  which  is  a  physically  meaningful  and  unique  state. 

5.3.1  Grid  Independence  Study 

The  computational  domain  in  this  study  is  chosen  as  2000x400/^m.  Three 

different  sets  of  grids  have  been  tested  to  determine  an  appropriate  computational  mesh. 

The  droplet  lifetimes  based  on  these  grids  are  compared,  as  shown  in  Table  5.1,  where 

the  incoming  hydrogen  velocity  is  20m/s  at  pressure  of  lOOatm.  Results  based  on  the  last 

two  sets  of  grids  are  very  close,  with  a  relative  error  about  4%.  Therefore,  the  second  set 

of  grids,  201x71 ,  is  considered  sufficient  for  this  numerical  study. 

T^l^ST.  Droplet  Lifetimes  at  Different  Numerical  Grids 

169x47  ~2oT^71  253x99' 

Droplet  Lifetimes  )  59  70  73 

Relative  Error  N/A  18.6%  4.3% 

5.3.2  Flow  Patterns 

The  flow  fields  are  illustrated  in  Figs.  5. 3-5. 6,  where  the  reference  frame  is 
moving  with  the  droplet  at  the  mass-weighted  droplet  velocity.  Figs.  5. 3-5. 5  present 
results  corresponding  to  a  pressure  of  100  atm,  and  Figs.  5. 5-5. 6  represent  flows  with  a 
same  Reynolds  number  (Re  =  120),  but  different  pressures.  In  all  these  figures. 
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streamlines  are  presented  with  temperature  contours,  which  is  used  to  indicate  the 
location  of  the  droplet.  It  should  be  emphasized  that  in  the  unsteady  fluid  flows,  the 
streamlines  could  only  represent  the  instant  velocity  directions.  Figure  5.3  illustrates  the 
transient  variations  of  the  flow  fields  at  an  incoming  velocity  of  2.5  m/s  and  a  pressure  of 
lOOatm  with  the  corresponding  Reynolds  number  (Re)  about  15.  In  this  low  Reynolds 
number  case,  there  is  no  recirculating  wake  behind  the  droplet,  and  the  incoming  flow  is 
redirected  and  passes  the  droplet  smoothly.  Because  of  the  droplet  vaporization,  there  is  a 
stagnant  point  in  front  of  the  droplet,  which  exists  in  the  entire  droplet  lifetime.  With  the 
incoming  velocity  increasing  to  5.0  m/s  and  Reynolds  number  to  30  as  shown  in  Fig.  5.4, 
the  flow  patterns  change  dramatically.  A  recirculating  eddy  appears,  and  its  center  moves 
downstream  during  the  vaporization  process.  When  compared  with  Fig  5.3,  the  stagnation 
point  in  front  of  the  droplet  moves  closer  to  the  droplet  surface  because  of  stronger 
incoming  aerodynamic  force.  Figure  5.5  corresponding  to  a  velocity  of  20  m/s  at  100  atm 
( Re  =  120 ),  where  a  recirculating  eddy  behind  the  droplet  is  visible  just  shortly  after  the 
droplet  is  injected.  The  eddy  is  moving  away  from  the  droplet,  and  it  is  growing  bigger. 
Results  in  Figs.  5. 3-5. 5  clearly  demonstrate  that  the  Reynolds  number  presents  strong 
effect  on  the  supercritical  fluid  flows  at  a  same  pressure  of  100  atm. 

In  Fig.  5.6,  pressure  is  increased  to  400  atm.  With  the  decrease  of  the  incoming 
velocity  to  5  m/s,  the  corresponding  Reynolds  number  is  about  the  same  as  that  in  Fig. 
5.5  (Re  =  120).  At  different  incoming  velocities,  the  same  Reynolds  number  is 
maintained  mainly  by  the  increase  of  the  hydrogen  density  with  pressure  as  shown  in 
Table  5.2a,  where  the  pressure  effects  on  thermodynamic  properties  are  presented.  For 
example,  the  hydrogen  density  is  2.4kg/m^  at  lOOatm  and  lOOOK,  but  it  becomes  slightly 
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above  9kgW  at  the  same  temperature  but  400atm.  At  the  same  Reynolds  number  but 
different  pressures,  the  streamlines  in  Fig.  5.6  look  similar  to  those  in  Fig.  5.5,  but  the 
recirculating  eddy  becomes  smaller.  Results  in  Figs  5. 5-5. 6  demonstrate  the  significant 
effect  of  pressure  on  flow  fields. 

In  fluid  mechanics  literature,  there  are  two  distinct  explanations  about  the 
formation  of  recirculating  wakes  behind  bluff  bodies  (Dandy  and  Leal  1986,  Leal  1989): 
vorticity  accumulation  and  boundary-layer  separation.  Leal  (1989)  argues  that  the 
recirculating  wake  behind  a  bluff  body  with  a  smooth  slip  surface  is  created  by  vorticity 
accumulation.  The  flow  detachment  occurs  when  the  maximum  surface  vorticity  reaches 
a  critical  value  about  5.25  (which  should  be  10.5  corresponding  to  the  definition  of  the 
dimensionless  vorticity  in  this  work,  co'  =  (O-dQjuQ).  Figure  5.7  presents  the  vorticity 

distributions  corresponding  to  the  four  cases  discussed  earlier.  The  maximum  vorticities 
in  all  these  cases  are  far  less  than  10.5,  especially  for  the  second  case  at  Re  =  30 .  The 
formation  of  the  recirculating  wake  in  this  supercritical  droplet  vaporization  case  is 
caused  by  flow  separation,  which  is  enhanced  by  droplet  deformation  and  surface 
blowing.  The  streamlines  in  Figs.  5. 3-5. 6  clearly  illustrate  that  the  flow,  resulting  from 
vaporization,  divides  the  flow  field  into  the  inner  and  outer  regions,  and  the  surface 
blowing  results  in  flow  separation  from  the  droplet  surface.  When  the  incoming  flow  is 
strong,  a  recirculating  eddy  appears.  The  further  stretching  of  the  droplet  perpendicular  to 
the  flow  direction  enhances  the  flow  separation,  and  helps  the  recirculating  eddy  grow 
bigger. 

In  Figs.  5. 3-5. 6,  the  flow  fields  inside  the  droplet  regions  clearly  illustrate  the 
trend  of  droplet  deformation,  which  will  be  discussed  shortly. 
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5.3.3  Temperature  Contours  and  Droplet  Deformations 

Under  supercritical  conditions,  an  oxygen  droplet  becomes  a  dense  fluid,  and  the 
droplet  surface  is  defined  as  an  isothermal  line  at  the  oxygen  critical  temperature. 
Therefore,  study  of  the  temperature  fields  will  reveal  the  droplet  deformation  and  the 
possible  secondary  break-up  phenomena,  which  are  closely  related  to  the  fuel  vapor 
distributions  and  mixing  in  the  combustion  devices.  Many  research  works  were 
conducted  about  droplet  deformations  and  secondary  break-ups  at  low  pressures  without 
heat  transfer  (Hinze  1955,  Range  and  Nicholls  1969,  Wierzba  and  Takayama  1988,  Pilch 
and  Erdman  1987,  Hsiang  and  Faeth  1992).  The  studies  have  determined  that  droplet 
deformation  and  break-up  patterns  in  a  forced  convection  flow  are  affected  by  two 
dimensionless  parameters:  the  Weber  number  and  the  Ohnesorge  number.  The  Weber 

number.  We  =  p^d^UQ  /a ,  represents  the  ratio  of  aerodynamic  force  to  the  surface 

tension,  and  the  Ohnesorge  number,  Oh  =  ,  is  the  ratio  of  viscous  force  to 

surface  tension.  Here  o  is  surface  tension,  dq  is  the  initial  droplet  diameter,  and  the 

subscripts  and  d  refer  to  hydrogen  and  oxygen,  respectively.  Because  the  oxygen 
droplet  is  assumed  to  reach  critical  mixing  state  immediately  after  being  injected  into  the 
combustion  chamber,  there  is  no  surface  tension  in  this  case.  One  reasonable  way  is  to 
define  another  dimensionless  parameter  in  this  supercritical  droplet  deformation  case  to 
represent  the  ratio  of  aerodynamic  to  viscous  forces. 

p^UQdp  p, 

Pd  I  -yj P^dQG  p^  ^ p^  Pd  Pd  y  Ps 


(5.1) 
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This  dimensionless  parameter  is  related  to  Reynolds  number,  hydrogen  and 
oxygen  viscosity  ratio,  and  square  root  of  density  ratio,  which  are  the  crucial  factors 
affecting  droplet  deformations  at  supercritical  conditions. 

To  illustrate  the  oxygen  droplet  deformation  and  the  effects  of  the  dimensionless 

parameter  Re-^  j-^,  the  temperature  contours  under  different  conditions  are 
/^d\Ps 

presented  in  Figs.  5. 8-5.1 1,  where  the  transient  developments  are  illustrated. 

Figures  5.8-5.10,  where  the  Reynolds  numbers  are  the  same  as  those  in  Figs.  5.3- 
5.5,  illustrate  the  Reynolds  number  effects  at  lOOatm.  Figure  5.8  presents  the  temperature 
contours  in  an  incoming  hydrogen  flow  at  a  low  velocity,  u=2.5m/s  (Re=15).  The  droplet 
slightly  extends  in  the  cross-flow  direction  and  becomes  an  oblate.  When  the  incoming 
velocity  increases  to  5m/s  (Re=30),  the  droplet  deformation  becomes  stronger  in  Fig.  5.9. 
The  droplet  first  becomes  an  oblate.  With  the  back  further  flattened  up,  the  droplet  begins 
to  look  like  a  ‘mushroom’.  The  stripping  effect  at  the  tip  of  the  droplet  appears  during  the 
vaporization  process.  As  the  Reynolds  number  increases  to  about  120  in  Fig.  5.10,  the 
strong  stripping  effect  is  clearly  illustrated,  and  a  sharp  tip  is  stretched  out.  The  droplet  is 
further  extended  in  the  cross-flow  direction  and  becomes  slender.  Afterwards,  it  is  easily 
bent  backward,  and  remains  a  crescent  shape  near  the  end  of  the  vaporization  process. 

The  pressure  effect  on  the  droplet  deformations  is  studied  in  Figs.  5.10-5.11, 
corresponding  to  a  same  Reynolds  number  ( Re  «  120 )  but  with  two  different  pressures  of 
lOOatm  and  400atm,  respectively.  As  expressed  in  Eq.  (5.1),  pressure  affects  the  droplet 
deformation  by  changing  ratios  of  viscosity  and  the  square  root  of  density.  As  shown  in 

Table  5.2b,  the  values  of  EA.  decrease  with  the  increase  of  pressure.  This  indicates 

Pd  V  Ps 
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the  droplet  deformation  will  become  weaker  at  a  higher  pressure.  The  results  in  Figs. 
5.10-5.11  confirm  this  conclusion.  With  the  increase  in  pressure,  the  stripping  effect  at 
the  tip  of  the  droplet,  as  well  as  the  extension  in  the  cross-flow  direction,  obviously 
weaken.  Fig.  5.11  shows  no  stripping  effect. 

Table  5.2a.  Pressure  Effects  on  Thermophysical  Properties  of  l^drogen  and  Oxygen 


Pressures 

(atm) 

100 

Hydrogen  at  lOOOK 
Density  Viscosity 

Oxygen  at  lOOK 
Density  Viscosity 

2.4 

( N  ■  s 1 ) 
1.968e-05 

(kglm^) 

1118.2 

(N-s/m^) 

1.721e-04 

200 

4.7 

1.972e-05 

1 144.9 

1.884e-04 

400 

9.0 

1.982e-05 

1184.3 

2.187e-04 

Table  5.2b.  Pressure  effects  on  the  Combined  Parameters 


Pressure  (atm) 

Pi.  Ip^ 

Pd 

2 

PdiPs 

Ps 

PR 

100 

2.47 

206.7 

200 

1.63 

184.2 

400 

1.04 

161.2 

The  analysis  above  demonstrates  that  the  ratio  of  aerodynamic  and  viscous  forces, 

as  defined  by  the  dimensionless  parameter  !  Oh  (=Re-^  l-^),  is  the  major 

^d\Ps 

factor  determining  the  droplet  deformation  under  supercritical  conditions.  At  a  given 
pressure,  the  variations  of  incoming  velocity  affect  the  droplet  deformation  by  changing 
the  Reynolds  number.  At  a  given  Reynolds  number,  pressure  will  affect  the  droplet 

deformation  by  changing  the  value  of  — 

Pd 


Pd. 

Ps 
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From  these  figures,  what  is  also  clear  is  that  the  droplet  shape  is  closely  related  to 
the  formation  of  the  recirculating  wake  behind  the  droplet,  as  shown  in  Figs.  5. 3-5. 6. 
With  the  forward  flow  close  to  the  axisymmetric  axis,  the  temperature  gradient  becomes 
very  steep  in  that  region. 

5.3.4  Mass  Fraction  Distributions 

The  mass  fraction  distributions  of  the  oxygen  species  are  shown  in  Fig.  5.12, 
where  the  same  cases  with  both  the  Reynolds  number  and  pressure  effects,  as  shown  in 
Figs.  5. 3-5. 6,  are  presented.  In  this  figure,  the  vaporized  oxygen  species  are  strongly 
advected  downstream  and  mixed  with  hydrogen  by  the  incoming  flow.  The  forward 
diffusion  is  limited  within  a  very  thin  boundary  layer,  and  the  diffusion  in  the  cross-flow 
direction  is  also  restricted  by  the  flow. 

Comparing  to  the  figures  of  temperature  contours,  the  core  region  of  the  oxygen 
species  deforms  in  the  same  patterns  with  the  strong  stripping  at  the  tip  and  extension  in 

the  cross-flow  direction  at  larger  values  of  ! Oh .  In  fact,  the  conclusions  can  also 

be  applied  to  the  droplet  region  defined  by  density  contours.  This  similarity  demonstrates 
that,  no  matter  how  we  define  the  droplet  surface,  the  droplet  deformation  is  mainly 

determined  by  the  dimensionless  parameter  We^^^/Oh  in  the  supercritical  droplet 
vaporization  processes. 

It  should  be  mentioned  that,  there  is  no  secondary  droplet  break-up  observed  in 
this  study  because  of  the  lack  of  surface  tension.  This  is  in  consistency  with  the  low- 
pressure  droplet  break-up  results  (Pilch  and  Erdman  1987,  Hsiang  and  Faeth  1992), 
which  indicate  that  break-up  seems  impossible  when  the  Ohnesorge  number 
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approaches  infinity.  This  is  exactly  the  case  here  with  no  surface 


tension. 

5.3.5  Droplet  Vaporization  Times 

In  spray  combustion  models,  the  droplet  vaporization  rate  and  droplet  lifetime 
must  be  determined.  The  droplet  vaporization  rates  at  supercritical  conditions  are 
presented  in  Figs.  5.13-5.15,  where  both  the  pressure  and  Reynolds  number  effects  are 
illustrated.  In  these  figures,  the  definition  of  the  droplet  mass  ratio  is 


Droplet  Mass  Ratio  = 


The  Remaining  Droplet  Mass  (m) 
The  Initial  Droplet  Mass  (mg) 


(5.2) 


Figure  5.13  shows  the  results  at  a  pressure  of  lOOatm  with  velocities  ranging  from 
2.5-20m/s,  which  correspond  to  the  Reynolds  numbers  (Re)  from  15-120.  At  the 


beginning  of  all  the  droplet  vaporization  processes,  there  are  rapid  droplet  vaporization 
periods,  because  droplets  (at  lOOK)  are  suddenly  encountering  the  high-temperature 
environments  (at  lOOOK).  After  very  short  times,  vaporization  rates  slow  down,  and 


droplets  at  different  Reynolds  numbers  begin  to  evaporate  at  different  rates.  The  droplet 
facing  a  higher  incoming  velocity  (Reynolds  number)  begins  the  transition  earlier 
because  of  the  stronger  convection  effect.  Figure  5.14  represents  the  same  situations  at 
400atm,  where  the  differences  among  the  transition  times  are  smaller,  because  of  the 


relative  higher  Reynolds  numbers  (60-230)  compared  with  the  cases  at  the  same 
velocities  in  Fig.  5.13.  The  results  in  these  two  figures  clearly  demonstrate  the  strong 
effects  of  Reynolds  numbers  on  the  droplet  vaporization  rates  and  droplet  lifetimes. 
Droplets  evaporate  much  faster  and  result  in  much  shorter  lifetimes  at  higher  Reynolds 


numbers. 


Time  (|lis) 

Fig.  5.13  Effects  of  incoming  velocities  (Reynolds  numbers)  on  the  droplet  vaporization 
rates  at  a  pressure  of  lOOatm. 


Time  (|lls) 

Fig.  5.14  Effects  of  incoming  velocities  (Reynolds  numbers)  on  the  droplet  vaporization 
rates  at  a  pressure  of  400atm. 
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At  a  given  Reynolds  number,  the  transition  to  convective  vaporization  starts 
earlier,  and  the  droplet  lifetime  is  shorter  at  a  lower  pressure  case.  Pressure  affects  the 
droplet  lifetimes  by  changing  the  Prandlt  and  Schmidt  numbers  of  the  oxygen-hydrogen 
mixture,  as  shown  in  Table  5.3,  where  the  properties  of  the  oxygen  and  hydrogen  mixture 
at  a  temperature  of  200K  and  equal  molar  fractions  are  presented  as  a  representative  case 
at  the  near-droplet  regions.  The  Schmidt  number  Sc  increases  significantly  with 
pressures,  but  the  Prandlt  number  Pr  decreases  slightly.  Therefore,  the  droplet  lifetimes 
are  proportional  to  the  Schmidt  number  Sc),  and  inversely  proportional  to  the  Prandlt 
number  ( ~  l/Pr ). 


Table  5.3.  Properties  of  Oxygen  and  Hydrogen  Mixture 
_ at  Equal  Molar  Fractions  and  200K  _ 


p  (atm) 

P  , 

Cp 

K 

p*10'' 

Dij=*=10^ 

Pr 

Sc 

Le 

(kg/m^) 

(J/kg-K) 

(W/m-K) 

(N-s/m^) 

(m^/s) 

100 

10T76 

2006.04 

0.054 

0.1489 

0.3877 

0.555 

0.370 

0.667 

200 

195.64 

2192.15 

0.070 

0.1728 

0.1984 

0.544 

0.445 

0.818 

400 

325.81 

2270.10 

0.098 

0.2250 

0.0992 

0.523 

0.696 

1.331 

Table  5.4.  Droplet JLifetimes  in  a  Quiescent  Environment 
Pressure  (atm)  100  200  400 

^  (Ps)  "  . 
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Results  in  Figs.  5.13-5.14  demonstrate  that  droplet  lifetime  is  a  function  of 
Reynolds  number  (Re)  and  pressure  at  a  supercritical  condition.  Before  attempting  any 
effort  to  find  the  possible  relationship,  the  variables  of  droplet  lifetime  ( t )  and  pressure 
are  first  non-dimensionalized.  An  apparent  choice  of  the  reference  pressure  is  the  critical 
pressure  of  oxygen  (Pc,02  =49.74ann).  The  reference  time  ((q)  is  chosen  to  be  the 

droplet  lifetime,  with  an  isolated  droplet  vaporizing  in  a  quiescent  environment  at  a  given 
pressure,  as  presented  in  Table  5.4.  Droplet  lifetimes  in  Table  5.4  are  calculated  from  the 
method  presented  in  Chapter  2. 

Utilizing  the  dimensionless  parameters,  the  following  relationship  is  what  we  are 
searching  for: 

-  =  /(Re,p^)  (5.3) 

^0 

where  is  the  dimensionless  pressure,  pj^  =  plPc,02  • 

At  given  pressures,  the  following  formula  is  a  natural  choice  following  the 
popular  Ranz  and  Marshall  correlation: 

= - 1 -  (5.4) 

^0  l-t-cRe° 

Equation  (5.4)  naturally  set  the  time  ratio  to  1  with  no  convection  ( Re  =  0 ). 

Data  Analysis  reveals  that  the  coefficient  a  is  approximately  a  constant  of  1.26  at 
different  pressures.  This  leads  to  the  conclusion  that  the  coefficient  c  must  be  a  funetion 
of  pressure.  The  following  relation  is  proposed: 

c  =  bpi 


(5.5) 
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Equation  (5.5)  does  not  change  the  unique  characteristics  of  Eq.  (5.4) 

(— =  1,  at  Re  =  0). 

^0 


Based  on  the  available  numerical  data,  the  following  coefficients  are  obtained: 
a  =  1.26, =  -1.58,6  =  0.1549  (5.6) 


Fig.  5.15  Correlation  of  the  dimensionless  droplet  lifetime  as  a  function  of  the  Reynolds 
number  and  pressure 
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Finally,  a  relationship  for  calculating  droplet  lifetimes  becomes 


t 

^0 


_ 1 _ 

1  + 0.1549 Re^-^^ 


(15  <  Re  <  240,2  <  R/?  <  8) 


(5.7) 


The  results  from  both  numerical  computations  and  Eq.  (5.7)  are  compared  in  Fig. 
5.15,  where  excellent  agreement  is  reached. 

5.3.6  Droplet  Dynamics 

In  this  section,  droplet  drag  coefficient  during  vaporization  process  is  studied.  In  a 
convective  environment,  a  droplet  constantly  deforms.  In  this  study,  the  droplet  location 
is  defined  as  its  center  of  mass.  Because  of  the  axisymmetry,  droplet  locations  always  lie 
in  the  axis  and  are  moving  downstream.  The  droplet  velocity  is  calculated  by  the  transient 
variations  of  droplet  locations  based  on  the  same  method  used  by  the  experimentalists 
(Yuen  and  Chen  1976).  Figures.  5.16-5.17  present  the  temporal  variations  of  droplet 
velocities,  and  the  effects  of  both  Reynolds  number  and  pressure  are  clearly  illustrated. 
Velocities  of  the  droplets  increase  with  time  in  both  figures.  At  a  given  pressure  of 
lOOatm  in  Fig.  5.16,  a  droplet  moves  faster  with  a  higher  incoming  velocity  (higher 
Reynolds  number).  At  a  given  Reynolds  number  of  60  in  Fig.  5.17,  a  droplet  moves  faster 
at  a  lower  pressure  because  of  the  faster  incoming  velocity.  Flowever,  there  are 
oscillations  in  both  figures,  especially  at  the  beginning  and  near  the  end  of  the 
vaporization  processes.  These  oscillations  are  caused  by  both  droplet  deformation  and 


surface  blowing. 
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Fig.  5.16  Transient  variations  of  droplet  velocities  and  the  linear  relationships  at  a 
pressure  of  lOOatm. 


Fig.  5.17  Transient  variations  of  droplet  velocities  and  the  linear  relationships  at  a 
Reynolds  number  of  60. 
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Due  to  the  sudden  injection  of  a  droplet  into  a  uniform  flow  field,  the  time 
required  for  the  flow  field  to  relax  can  be  estimated  as  =  r/o /mq-  Careful 

inspection  of  the  transient  variations  of  droplet  velocities  indicates  simple  linear  relations, 
except  within  the  relaxation  times  when  the  nonlinear  effects  are  trivial  because  of  the 
relative  slow  velocities  and  short  times.  Therefore,  the  transient  variations  of  droplet 
velocities  are  correlated  as 

V  =  bt  (5-8) 

It  should  be  kept  in  mind  that  the  coefficient  b  has  a  dimension  of  m/s  .  The 
following  relationship  of  the  velocity  coefficient  b ,  which  is  a  function  of  Reynolds 
number  and  pressure,  is  assumed: 

b  =  cRt^pi  (5-9) 

In  Eq.  (5.9),  when  the  Reynolds  number  is  set  to  0,  a  zero  value  of  coefficient  b 
will  naturally  lead  to  a  stationary  droplet,  as  physically  required. 

All  the  parameters  in  Eq.  (5.9)  can  be  easily  determined  by  simple  linear 
regression.  Analyzing  the  available  data  at  given  pressures  indicates  that  the  coefficient 
a  is  a  constant  of  1.76.  At  given  Reynolds  numbers,  the  coefficient  d  also  approaches  a 
constant.  Finally,  with  the  known  coefficients  a  and  ,  the  coefficient  c  can  be  easily 
found  out.  These  coefficients  are 

a  =  1.76,r/  =  -1.164,c  =  12  (5.10) 

where  the  dimension  of  parameter  c  is  m/sl  Substituting  these  numbers  into  Eq.  (5.9), 
the  relationship  for  the  velocity  coefficient  b  is  expressed  as 

b  =  12  Re^-^^  (15  <  Re  <  240,2  <Pr<^) 


(5.11) 
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Figures  5.18  shows  a  very  good  agreement  between  the  numerical  results  and 
those  from  Eq.  (5.11). 

The  correlation  in  Eq.  (5.11)  is  employed  to  calculate  the  transient  variations  of 
droplet  locations,  as  illustrated  in  Figs.  5.19-5.20.  The  droplet  location  S  is  calculated  as 

S=-bt^  (4.12) 

2 


Fig.  5.18  Correlation  of  the  velocity  coefficient  b  as  a  function  of  the  Reynolds  number 
and  pressure. 
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Fig.  5.19  Transient  variations  of  droplet  locations  from  both  correlation  and  numerical 
computations  at  a  pressure  of  lOOatm. 


Fig.  5.20  Transient  variations  of  droplet  locations  from  both  correlation  and  numerical 
computations  at  a  Reynolds  number  of  60. 
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Figures  5.19-5.20  represent  two  cases  at  /?  =  lOOafm  and  Re  =  60,  respectively. 
In  these  figures,  droplet  locations  calculated  from  Eq.  (5.12)  are  compared  with  those 
directly  from  numerical  computations.  Good  agreements  are  clearly  demonstrated. 

Based  on  Eq.  (5.8),  it  is  easy  to  find  the  droplet  acceleration  , 

a,=b  (5.13) 

The  drag  force  acting  on  the  droplet  can  be  approximately  expressed  as 
FD=ma^  (5.14) 

The  drag  coefficient  is  defined  as 


Cd  - 


D 


lllfyulnrl 


(5.15) 


Because  of  the  droplet  deformation,  the  droplet  radius  r„■^  is  defined  as  the  radius 


of  a  mass-equivalent  droplet  sphere  with  a  uniform  density  at  its  initial  value. 


f 


'm 


3m 

ApdTt 


1/3 


(5.16) 


V  '  “  y 

where  m  is  the  droplet  mass  at  any  instant. 

Substituting  Eqs.  (5.14),  (5.16)  and  (5.11)  into  the  drag  coefficient  expression 
(5.15),  the  following  relationship  can  be  established; 


m 


Pd 

Ps 


Vr 


0^0.24  ,1.164 
Re  pji 


(15<Re<240,2<PR  <8) 


(5.17) 


where  is  defined  as  the  reference  kinematic  viscosity  of  hydrogen  at  lOOOK  and 


lOOatm. 
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Equation  (5.17)  clearly  indicates  that  the  drag  coefficient  always  decreases  with 
f  m 

time  ( ~  —  )  in  the  supercritical  vaporization  case.  This  trend  is  different  from  that 

of  the  low-pressure  droplet  vaporization  case  without  droplet  deformation  (Chiang  et  al. 
1992).  Results  from  Chiang  et  al.  (1992)  indicate  that  the  drag  coefficient  could  increase 
or  decrease  with  time  during  different  portions  of  the  droplet  lifetime.  At  the  end  of  the 
vaporization  process,  it  generally  increases  with  time  as  a  result  of  the  reduction  in  the 
Reynolds  number,  which  is  defined  using  the  relative  velocity  between  the  incoming 
stream  and  the  droplet.  In  the  supercritical  droplet  vaporization  case,  the  moving 
velocities  of  the  droplets  are  relatively  small  compared  with  the  incoming  velocity,  as 
shown  in  Figs.  5.16-5.17.  Therefore,  the  Reynolds  number  reduction  is  not  strong  enough 
to  increase  the  drag  coefficient. 

At  a  given  Reynolds  number,  pressure  influences  the  drag  coefficient  by  changing 


the  product  of 

Ps\^s] 


as  illustrated  in  Table  5.2b.  Results  in  Table  5.2b 


clearly  demonstrate  that  the  drag  coefficient  has  a  slight  decrease  with  the  increase  of 
pressure  at  a  given  Reynolds  number.  It  should  be  kept  in  mind  that  pressure  could  also 
affect  the  Reynolds  number  in  this  supercritical  droplet  vaporization  case. 


5.4  Summary  of  Results 

In  this  paper,  an  oxygen  droplet  vaporizing  in  convective  hydrogen  environments 
at  supercritical  conditions  are  numerically  investigated.  A  consistent  and  efficient 
numerical  algorithm  is  utilized,  with  full  account  of  real  fluid  behavior.  Thermodynamic 
relationships  are  derived  using  fundamental  theories  and  can  be  calculated  based  on  any 
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equation  of  state.  A  preconditioning  numerical  scheme  is  employed  to  handle  small- 
Mach-number  fluid  flows.  The  details  of  the  numerical  and  thermodynamic  treatments 
can  be  found  in  Chapter  4.  The  transport  properties  are  estimated  by  the  state-of-the-art 
techniques  with  proper  consideration  of  high-pressure  effects. 

An  oxygen  droplet  is  assumed  to  reach  the  critical  mixing  state  immediately  after 
being  injected  into  the  combustion  chamber.  Therefore,  this  is  a  problem  of  supercritical 
fluid  flow  with  strong  heat  and  mass  transfer.  An  isothermal  line  at  the  critical 
temperature  of  oxygen  is  defined  as  the  droplet  surface.  The  transient  developments  of 
temperature,  oxygen  mass  fraction,  and  flow  fields  are  presented.  The  effects  of  incoming 
velocity  (or  Reynolds  number)  and  pressure  on  droplet  deformations  are  explored  in 
detail.  The  major  factor  determining  the  droplet  deformations  under  supercritical 

conditions  is  a  dimensionless  parameter  We^^^/oh  (=Re-^  1-^^),  which  represents 

V  Ps 

the  ratio  of  aerodynamic  and  viscous  forces.  Both  Reynolds  number  and  pressure  perform 
their  effects  by  changing  this  parameter. 

The  droplet  vaporization  rates  and  droplet  lifetimes  are  studied.  Droplet  lifetimes 
are  affected  strongly  by  the  incoming  hydrogen  velocity  (or  Reynolds  number)  and 
pressure.  At  a  given  pressure,  a  droplet  evaporates  more  rapidly  at  a  higher  Reynolds 
number.  The  droplet  lifetime  is  significantly  shortened  at  a  lower  pressure  and  a  given 
Reynolds  number.  A  relationship  is  established  to  calculate  droplet  lifetimes  accurately. 

The  transient  variations  of  droplet  velocities  and  drag  coefficients  during  droplet 
vaporization  processes  are  illustrated.  A  linear  relationship  is  reasonably  assumed  for 
droplet  velocity,  and  the  linear  velocity  coefficient  b  is  found  to  be  proportional  to 
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p'r'^^  ■  Droplet  positions  calculated  from  this  linear  relationship  agree  well  with 
those  directly  from  numerical  computations.  The  drag  coefficient  always  decreases  with 


time  ( - 


—  )  as  a  result  of  the  weak  Reynolds  number  reduction  in  this  supercritical 


vaporization  case.  At  a  given  Reynolds  number,  pressure  influences  the  drag  coefficients 


by  changing  the  product  of 


Pd_ 

Ps 


and  the  drag  coefficient  decreases 


slightly  with  the  increase  of  pressure. 

Results  in  this  study  give  a  deeper  understanding  of  the  supercritical  droplet 
vaporization  process.  The  correlations  regarding  droplet  lifetime  and  drag  coefficient  can 
be  directly  applied  to  the  numerical  modeling  of  supercritical  spray  combustion. 
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Fig.  5.3  Transient  developments  of  flow  fields  at  u=2.5m/s,  p=100atm  (Re  =  15),  at 
instants  of  t=20, 100,  200, 400ps,  respectively. 
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Fig.  5.4  Transient  developments  of  flow  fields  at  u=5.0m/s,  p=100atm  (Re  =  30) ,  at 
instants  of  t=12,  100,  150,  250ps,  respectively. 
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Fig.  5.5  Transient  developments  of  flow  fields  at  u=20.0m/s,  p=100atm  (Re  =  120),  at 
instants  of  t=4,  20,  40,  60ps,  respectively. 


Fig.  5.8  Transient  developments  of  temperatu 
(Re  =  15) ,  at  instants  of  20,  100,  200, 400ps,  respi 


Fig.  5.9  Transient  developments  of  temperature  co; 
(Re  =  30) ,  at  instants  of  12, 100,  150, 250ps,  respective! 


Fig.  5.10  Transient  developments  of  temperati 
(Re  =  120) ,  at  instants  of  4, 20, 40,  60ps,  respect: 
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Fig.  5.12  Distributions  of  oxygen  mass  fractions,  (a)  u=2.5ni/s,  p=100at 
u=5.0m/s,  p=100atm,  t=150)is,  (c)  u=20.0m/s,  p=100atm,  t=40ps, 
p=400atm,  t=60ps 


Chapter  6 


SUPERCRITICAL  VAPORIZATION  OF  TWO-DROPLETS-IN- 
TANDEM  IN  CONVECTIVE  ENVIRONMENTS 

6.1  Problem  Description 

The  effects  of  droplet  interactions  on  droplet  lifetimes  at  high  pressures  have  been 
studied  in  Chapter  3,  where  the  droplet  cluster  behavior  in  a  quiescent  environment  is 
extensively  examined.  However,  convection  effect  could  easily  change  the  results.  Since 
it  is  still  not  a  reality  for  current  computational  capacity  to  investigate  numerically  a 
group  of  droplet  vaporizing  in  a  convective  environment,  oxygen  droplet  interactions  in 
convective  hydrogen  environments  are  studied  in  this  chapter  by  setting  up  a  simple 
configuration  of  two  droplets  moving  in  tandem.  The  goal  is  to  identify  the  detailed 
vaporization  characteristics  with  two  neighboring  droplets  interacting  with  each  other, 
and  the  study  is  focused  on  the  influences  from  the  immediate  neighboring  droplet  on 
droplet  lifetimes  and  dynamics,  such  as  droplet  collision/separation  phenomenon.  The 
configuration  is  sketched  in  Fig.  6.1,  where  the  initial  droplet  spacing  is  prescribed  (H). 
The  ratio  of  the  initial  droplet  spacing  and  droplet  radius,  H/R,  is  defined  to  indicate  the 
intensity  of  droplet  interactions. 


Fig.6.1.  Two-droplets-in-tandem  in  a  convective  environment. 
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The  following  assumptions  are  taken  in  this  study: 

1)  The  flow  is  laminar  and  axisymmetric; 

2)  The  two  droplets  reach  the  critical  mixing  state  immediately  after  being 
injected  into  the  combustion  chamber; 

3)  The  droplet  boundary  is  defined  as  the  isothermal  line  at  the  critical 
temperature  of  oxygen  species  (154.6  K). 

6.2  Numerical  Treatment 

The  preconditioning  numerical  algorithm  with  a  unified  treatment  of  real  fluid 
behavior,  as  developed  in  Chapter  4,  is  utilized  to  solve  the  current  problem.  The 
methods  presented  in  Chapter  2  are  applied  to  evaluate  the  transport  properties,  with  full 
account  of  pressure  effect.  Non-uniform  grids  are  generated  in  the  same  manner  as  that  in 
the  single  droplet  case  so  that  no  extra  effort  is  needed  for  grid  independence  study.  Part 
of  the  grids  is  illustrated  in  Fig.  6.2,  where  meshes  around  the  two  droplets  are 
significantly  refined  for  numerical  accuracy. 


Fig.  6.2.  Part  of  the  Numerical  Grids  for  Two-Droplets-in-Tandem. 
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6.3  Results  and  Discussions 

Supercritical  vaporization  of  two-droplets-in-tandem  is  illustrated  in  Fig.  6.1, 
where  two  oxygen  droplets  are  suddenly  injected  into  a  uniform  hydrogen  flow  field. 
Initially,  the  two  droplets  are  at  the  same  size  of  50pm  in  diameter.  The  oxygen  droplet 
and  the  hydrogen  fluid  are  at  the  uniform  temperatures  of  100  K  and  lOOOK,  respectively. 
Pressures  range  from  100  to  400  atm.  As  for  the  single  droplet  case,  the  flow  is  sustained 
as  laminar  and  axisymmetric  by  limiting  the  Reynolds  Number  (Re)  to  be  less  than  240. 
The  Reynolds  number  is  defined  based  on  the  initial  droplet  diameter  and  the  inflow 
hydrogen  velocity,  density,  and  viscosity. 

6.3.1  Droplet  Deformations  and  Temperature  Fields 

As  in  the  single  droplet  vaporization  case,  the  droplet  boundary  is  defined  as  an 
isothermal  line  at  the  critical  temperature  of  the  oxygen  species  (154.6  K)  in  the 
numerical  investigation.  Therefore,  the  temperature  fields  will  illustrate  the  droplet 
deformation  characteristics. 

Figures  6. 3-6. 5  present  the  typical  droplet  deformation  features  at  a  pressure  of 
lOOatm,  an  incoming  velocity  of  lOm/s  (Re  =  60),  but  different  H/R  ratios  of  4,  8,  and 
12.  Figure  6.3  shows  the  transient  developments  of  temperature  fields  at  a  H/R  ratio  of  4. 
Initially,  the  two  droplets  evaporate  independently.  With  the  increase  of  time,  the 
convection  effect  begins  to  dominate  the  vaporization  process,  and  the  temperature  fields 
around  the  two  droplet  merges.  At  the  time  of  t  =  60  iJ.s ,  the  two  droplets  deform  in 
completely  different  patterns.  The  leading  droplet  confronting  stronger  incoming  flow 
behaves  more  like  an  isolated  droplet  and  deforms  perpendicular  to  the  flow  direction.  Its 
backside  becomes  flat,  while  the  front  side  retains  a  round  shape.  The  deformation  of  the 
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trailing  droplet,  due  to  much  weaker  incoming  fluid  flow,  is  dominated  by  diffusion.  The 
phenomenon  is  especially  profound  near  the  axisymmetric  axis,  where  the  temperature 
contours  significantly  diffuse  outward,  resulting  in  relatively  small  temperature  gradient. 
The  boundary  layer  theory  is  obviously  invalid  in  this  region.  Near  the  axisymmetric 
axis,  the  trailing  droplet  strongly  expands  along  the  flow  direction  and  presents  the 
characteristics  of  deformation  distinct  from  the  leading  droplet.  When  time  is  increased  to 
155  /us,  the  leading  droplet  moves  very  close  to  the  trailing  droplet,  but  there  is  no 
droplet  collision  based  on  temperature  contours.  At  this  instant,  the  leading  droplet  is 
clearly  breaking  up  at  the  axisymmetric  axis.  With  most  of  mass  moving  toward  the  edge 
of  the  droplet,  the  break-up  is  categorized  as  the  typical  forward  bag  break-up.  The  mass 
ring  resulted  from  the  bag  break-up  is  transported  downstream  by  moving  over  the 
surface  of  the  trailing  droplet  and  disappears  shortly  after  the  break-up.  Meantime,  there 
is  apparent  stretching  at  the  tip  of  the  trailing  droplet.  The  outside  temperature  fields 
behave  as  if  they  are  surrounding  a  large  isolated  droplet.  After  the  leading  droplet 
disappears,  the  remaining  trailing  droplet  begins  to  deform  like  an  isolated  droplet,  except 
that  the  outside  temperature  contours  experience  stronger  expansion  perpendicular  to  the 
flow  direction.  Because  the  trailing  droplet  now  has  to  confront  incoming  flow  directly, 
strong  stretching  of  the  droplet  is  clearly  illustrated. 

Droplet  deformations  at  a  H/R  ratio  of  8  present  features  distinct  from  that  at  a 
H/R  ratio  of  4.  At  the  time  of  t  =  145  ,  the  leading  droplet  is  forced  to  expand  strongly 

perpendicular  to  the  axisymmetric  axis  and  retains  a  slightly  convexed  disk  shape. 
Meantime,  the  expansion  of  the  trailing  droplet  along  the  flow  direction  is  significantly 
weakened  compared  with  that  at  a  H/R  ratio  of  4,  although  it  is  still  dominated  by 
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diffusion  near  the  axisymmetric  axis.  The  trailing  droplet  is  strongly  bent  downstream 
and  shows  a  skirt-like  shape.  Before  the  low-temperature  core  regions  of  the  two  droplets 
contact  each  other,  the  leading  droplet  has  already  been  completely  evaporated,  as 
illustrated  at  the  instant  of  t=190  /us .  At  the  same  time,  the  trailing  droplet  is  further 
stretched.  There  is  no  droplet  break-up  occurring  for  the  leading  droplet  during  its  entire 
vaporization  process.  By  the  end  of  the  vaporization  process,  only  one  highly  stretched 
trailing  droplet  is  remaining  in  the  flow  field. 

Droplet  deformations  at  a  H/R  ratio  of  12  are  presented  in  Fig.  6.5.  During  the 
entire  droplet  lifetimes,  the  two  droplets  both  behave  in  similar  patterns  as  an  isolated 
droplet,  with  the  droplet  expansions  perpendicular  to  the  axisymmetric  axis  and  crescent 
shapes  remaining  at  the  end  of  the  vaporization  process. 

Increasing  the  incoming  velocity  to  20  m/s  (Re  ==120),  the  two  droplets  still 
deform  in  the  same  features  as  those  illustrated  in  Figs.  6.3-6. 5.  At  a  higher  pressure, 
droplet  interactions  become  weaker.  When  pressure  increases  to  400  atm,  the  two 
droplets  begin  to  deform  like  two  isolated  droplets  at  a  H/R  ratio  of  8,  while  the  same 
situation  only  occurs  at  a  H/R  ratio  of  12  at  a  pressure  of  100  atm. 

6.3.2  Mass  Fraction  Fields 

The  typical  mass  fraction  fields  are  presented  in  Figs.  6. 6-6. 7,  which  correspond 
to  the  same  initial  conditions  as  those  in  Figs.  6. 3-6.4.  At  a  H/R  ratio  of  4,  the  mass 
fraction  contours  around  the  two  droplets  begin  to  merge  even  at  the  very  beginning  of 
the  vaporization  process,  as  shown  in  Fig.  6.6.  When  time  increases  to  60  fis ,  the 
deformations  of  the  two  droplets  illustrated  by  the  mass  fraction  contours  are  similar  to 
those  represented  by  temperature  contours  shown  in  Fig.  6.3.  At  this  instant,  the  leading 
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droplet  becomes  an  oblate  with  a  flat  back,  while  the  trailing  droplet  maintains  a  prolate 
shape  with  indentation  at  its  backside.  The  difference  between  the  fields  of  mass  fraction 
and  temperature  is  that  species  diffuse  much  faster  than  thermal  diffusion  does,  and  there 
is  consequently  stronger  downstream  transport  of  the  oxygen  species.  In  Fig.  6.6,  as  time 
further  increases  to  125  jUs ,  the  two  droplets  collided  according  to  the  mass  fraction 
field,  and  a  complete  mushroom-like  shape  is  produced.  This  feature  is  significantly 
different  from  that  in  the  corresponding  temperature  field.  In  addition,  there  is  no  forward 
bag  break-up  observed  for  the  leading  droplet  in  the  mass  fraction  field.  At  the  end  of  the 
vaporization  process,  the  two  droplets  completely  merge  together,  and  begin  to  behave 
like  an  isolated  droplet. 

Fig  6.7  presents  the  developments  of  mass  fraction  fields  at  a  H/R  ratio  of  8.  No 
droplet  interaction  is  illustrated  with  the  vaporization  time  up  to  16  /«  .  As  time  reaches 
60  jiis ,  the  oxygen  species  from  the  leading  droplet  is  connected  with  those  from  the 
trailing  droplet  due  to  strong  convective  transport.  However,  at  this  time,  the  effect  of 
interactions  on  both  droplets  is  trivial,  and  the  two  droplets  still  perform  similar 
deformations.  With  time  increase  to  190  pis,  the  dense  oxygen  regions  of  the  two 
droplets  merge  together,  and  the  effect  of  droplet  interactions  becomes  significant.  The 
leading  droplet  expands  strongly  perpendicular  to  the  flow  direction,  while  the  trailing 
droplet  is  mainly  stretched  along  the  flow  direction.  However,  there  still  remains  the 
oxygen-rich  region  for  the  leading  droplet,  in  contrast  to  its  disappearance  in  the 
temperature  field.  Finally,  the  two  droplets  collide  and  form  one  droplet,  which  is  slightly 
stretched  at  the  tip.  The  stretched  deformation  of  the  merged  droplet  is  different  from  that 
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of  an  isolated  droplet,  which  maintains  a  crescent  shape  at  the  end  of  the  vaporization 
process. 

The  effect  of  droplet  interactions  on  the  deformations  of  the  two  droplets  becomes 
extremely  weak  at  a  H/R  ratio  of  12.  Both  droplets  deform  like  an  isolated  droplet  and 
remain  crescent  shapes  at  the  end  of  the  vaporization  process.  Increase  of  the  incoming 
velocity  from  10  m/s  to  20  m/s  does  not  change  the  deformation  characteristics.  In 
contrast,  increase  of  pressure  from  100  atm  to  400  atm  will  weaken  the  droplet 
interactions  just  like  the  situations  demonstrated  in  temperature  fields.  Only  minor  effect 
of  droplet  interactions  on  both  droplets  is  illustrated  once  the  H/R  ratio  reaches  8. 

6.3.3  Flow  Field  Developments 

With  two  droplets  confronting  the  forced  incoming  flow  in  tandem,  the  flow 
fields  around  the  leading  and  trailing  droplets  are  extremely  different.  The  characteristics 
distinct  from  those  of  an  isolated  droplet  are  also  displayed. 

Figures  6.8-6.10  present  the  flow  fields  at  the  same  initial  pressure  and  incoming 
velocity  of  100  atm  and  20  m/s  (Re  =  120)  but  different  H/R  ratios.  A  stationary 
reference  frame  is  maintained  in  these  three  figures.  Streamlines  are  utilized  to  illustrate 
the  instant  flow  fields  with  temperature  fields  displayed  to  indicate  the  droplet  locations. 
Figure  6.8  illustrates  the  transient  flow  developments  at  a  H/R  ratio  of  4.  At  the 
beginning  of  the  vaporization  process,  the  flow  fields  around  the  two  droplets  are  similar 
in  the  fact  that  a  recirculating  eddy  is  developed  in  the  wake  of  each  droplet.  Because  the 
leading  droplet  faces  stronger  incoming  flow  than  the  trailing  droplet  does,  flow 
separation  around  it  occurs  earlier,  resulting  in  a  slightly  larger  recirculating  region  in  its 
wake.  At  the  time  of  t  =  40  ,  as  the  leading  droplet  retains  a  disk-like  shape,  flow 
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separation  is  further  enhanced,  resulting  in  a  larger  recirculating  eddy.  Since  the  trailing 
droplet  diffuses  significantly  in  the  upstream  direction,  it  pushes  the  first  recirculating 
eddy  slightly  upward.  The  second  recirculating  eddy  also  grows  during  the  vaporization 
process,  but  it  is  relatively  smaller  compared  with  the  first  one.  The  reason  is  that,  with 
the  trailing  droplet  deforming  mainly  along  the  flow  direction,  flow  separation  around  it 
is  weaker  than  that  around  the  leading  droplet.  When  the  leading  droplet  moves  further 
close  to  the  trailing  droplet  (at  t  =  65  ^ ),  the  first  recirculating  eddy  is  pushed  upward 
and  located  at  the  top  of  the  trailing  droplet.  Because  the  presence  of  the  trailing  droplet 
restricts  the  possible  recirculating  space  for  the  first  eddy,  its  size  becomes  smaller.  With 
the  first  eddy  pushed  on  top  of  the  trailing  droplet,  an  interesting  phenomenon  appears. 
Since  flow  is  now  redirected  by  the  first  eddy  and  separated  behind  it  instead  of  the 
trailing  droplet,  the  original  second  recirculating  eddy  disappears.  Instead,  another 
recirculating  eddy  is  created  behind  the  first  one.  These  two  eddies  locate  very  close  to 
each  other.  At  the  time  of  t  =  85  fxs ,  the  two  droplets  move  so  close  to  each  other  that 
their  effect  on  the  flow  fields  become  similar  to  that  of  an  isolated  droplet.  Meantime,  the 
first  eddy  and  the  new  eddy  behind  it  merge  together,  resulting  in  one  large  recirculating 
region. 

The  flow  fields  at  a  H/R  ratio  of  8  is  different  from  that  at  a  H/R  =  4,  as  shown  in 
Fig.  6.9.  Initially,  there  are  two  eddies  independently  developing  behind  the  two  droplets. 
With  the  increase  of  time,  the  two  eddies  both  grow  in  size.  However,  the  first  eddy  is 
much  larger  than  the  second  one  because  of  two  reasons.  The  first  reason  is  related  to  the 
deformation  of  the  leading  droplet,  which  enhances  flow  separation.  The  second  reason  is 
that,  since  the  two  droplets  form  a  close  region  between  them,  it  helps  entraining  vorticity 
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and  maintaining  a  strong  recirculating  flow  there.  Once  the  leading  droplet  moves  closer 
to  the  trailing  droplet,  the  first  recirculating  eddy  is  pushed  upwards  and  decreases  in 
size.  The  size  of  the  second  eddy  also  decreases.  This  is  partially  due  to  the  deformation 
of  the  second  droplet,  and  partially  caused  by  the  position  of  the  first  eddy,  which 
significantly  changes  the  flow  field.  Since  the  leading  droplet  completely  disappears  at 
the  time  of  t  =  110  ,  the  result  is  the  disappearance  of  the  first  eddy.  However,  the 

second  eddy  is  maintained  during  the  entire  vaporization  process. 

The  flow  fields  corresponding  to  a  H/R  ratio  of  12  is  illustrated  in  Fig.  6.10, 
where  the  interactions  between  the  two  droplets  are  extremely  weak.  The  first 
recirculating  eddy  remains  growing  in  the  entire  vaporization  process  because  of  the 
existence  of  the  closed  region  between  the  two  droplets  and  its  entraining  of  vorticity. 
The  second  recirculating  eddy  increases  in  size  at  the  beginning  of  the  process,  but  its 
size  decreases  gradually  at  the  later  stage,  because  the  deformation  of  the  trailing  droplet 
eases  flow  separation. 

To  demonstrate  the  effects  of  incoming  velocity  and  pressure  on  flow  field 
developments.  Fig.  6.11  shows  the  flow  fields  at  a  H/R  ratio  of  8  for  two  different  initial 
conditions.  Fig  6.11a  corresponds  to  the  flow  field  at  a  pressure  of  100  atm  and  an 
incoming  velocity  of  lOm/s  (Re  =  60),  while  Fig  6.11b  represents  the  flow  fields  at  a 
pressure  of  400  atm  and  an  incoming  velocity  of  u  =  5  m/s  ( Re  =  120 ).  The  flow  field  in 
Fig.  6.11a  indicates  that  the  second  recirculating  eddy  disappears  without  creation  of  a 
new  recirculating  region,  because  incoming  flow  is  not  sufficiently  strong  to  maintain  a 
second  eddy.  Figure  6.11b  corresponds  to  a  same  Reynolds  number  (Re  =  120)as  Fig. 
6.9  does,  but  its  flow  field  is  completely  different  from  the  former  one  in  the  fact  that 
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there  are  two  independent  recirculating  eddies  existing  at  the  end  of  the  vaporization 


process.  Comparison  of  the  flow  fields,  as  shown  in  both  Fig. 6. 9  and  6.11b,  indicates  that 
droplet  interaction  becomes  significantly  weaker  at  a  higher  pressure  of  400  atm.  This 


conclusion  is  consistent  with  those  drawn  from  both  the  temperature  and  mass  fraction 
fields. 


a) 


£ 

6 


X  (mm) 


910.5 
820.8 

731.2 

641.5 

551.9 

462.2 

372.6 

282.9 

193.3 
155.0 


b) 


B 

B 


X  (mm) 


M  910.9 
^  821.7 

—  732.6 

—  643.4 
554.3 
465.1 

—  376.0 
h-i  286.8 


Fig.  6.11  Effects  of  Incoming  Velocity  and  Pressure  on  the  Flow  Fields,  a)  p=100  atm, 
u=10  m/s,  t=145ps,  b)  p=400  atm,  u=5  m/s,  t=115  ps. 

6.3.4  Droplet  Lifetimes 

With  two  droplets  moving  in  tandem,  the  droplet  vaporization  rates  and  droplet 
lifetimes  for  the  leading  and  trailing  droplets  show  significantly  different  characteristics. 
Figures  6.12-6.13  illustrate  the  vaporization  rates  for  the  two  droplets  at  different  H/R 
ratios.  Figure  6.12  represents  the  cases  at  an  initial  pressure  and  velocity  of  100  atm  and 
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20  m/s  (Re  =  120),  respectively.  These  figures  clearly  demonstrate  that  the  droplet 
interactions  at  different  H/R  ratios  present  only  minor  effects  on  the  leading  droplets,  but 
effects  on  the  trailing  droplets  are  much  stronger.  The  trailing  droplet  evaporates  at  the 
same  rates  as  the  leading  droplet  only  at  the  beginning  of  the  vaporization  process.  Once 
the  two  droplets  start  to  interact  with  each  other,  the  vaporization  rates  of  the  trailing 
droplets  begin  to  deviate  significantly  from  those  of  the  leading  droplets.  The  exact 
deviation  times  depend  on  the  H/R  ratios.  At  a  H/R  ratio  of  4,  the  vaporization  rate  of  the 
trailing  droplet  slows  down  only  shortly  after  the  beginning  of  the  vaporization  process. 
During  the  times  around  90-100  jiis ,  there  is  another  sudden  change  of  the  vaporization 
rate.  During  this  second  transition  period,  the  leading  droplet  breaks  up.  The  mass  ring 
resulted  from  the  break-up  is  transported  over  the  surface  of  the  trailing  droplet  and 
causes  the  low-temperature  oxygen-rich  region  of  the  leading  droplet  to  merge  with  the 
trailing  droplet.  There  is  no  obvious  transition  period  existing  in  the  other  two  cases  at  the 
H/R  ratios  of  8  and  12,  since  the  low-temperature  regions  of  the  leading  droplets  have 
already  disappeared  before  they  could  contact  the  trailing  droplets. 

At  a  pressure  of  400  atm  and  an  incoming  velocity  of  5  m/s  (Re  =  120),  droplet 
interaction  effects  on  droplet  vaporization  rates  become  extremely  weak  as  demonstrated 
in  Fig.  6.13.  Significant  change  of  droplet  vaporization  rates  of  the  trailing  droplet  only 
occurs  when  the  two  droplets  are  very  close  to  each  other  at  a  H/R  ratio  of  4. 

To  further  clarify  the  effect  of  droplet  interactions  on  droplet  lifetimes,  the  droplet 
lifetimes  of  both  the  leading  and  trailing  droplets  are  compared  with  that  of  an  isolated 
droplet,  while  all  the  other  flow  conditions  are  maintained  as  the  same.  A  droplet  lifetime 


ratio  7  is  defined  as 
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Lifetime  of  an  isolated  droplet  , ,  ^ . 

T  = -  (0.1) 

Lifetime  of  the  leading/trailing  Droplet 

Figure  6.14-6.15  present  the  variations  of  droplet  lifetime  ratios  at  two  pressures 
of  100  and  400  atm  and  two  incoming  velocities  of  10  and  20m/s.  At  a  pressure  of  100 
atm,  there  is  only  slight  interaction  effect  on  the  droplet  lifetimes  of  the  leading  droplet  at 
a  H/R  ratio  of  4,  at  which  the  droplet  lifetime  ratios  i  are  around  0.8  for  both  incoming 
flow  velocities.  The  lifetimes  of  the  trailing  droplets  are  dramatically  elongated  due  to  the 
weak  convection  effect  caused  by  the  existence  of  the  leading  droplet.  The  droplet 
lifetime  ratios  are  about  the  same  (around  0.5)  for  the  trailing  droplets  at  both  H/R  ratios 
of  4  and  8,  although  their  initial  vaporization  rates  are  different.  At  a  H/R  ratio  of  4,  the 
trailing  droplet  initially  evaporates  at  a  slower  rate,  but  the  leading  droplet  in  this  case 
disappears  in  relatively  shorter  time  than  that  at  a  H/R  ratio  of  8.  After  the  disappearance 
of  the  leading  droplet,  the  trailing  droplet  begins  to  evaporate  like  an  isolated  droplet,  and 
its  vaporization  rate  consequently  becomes  faster  than  before. 

Figure  6.14  clearly  indicates  that  there  are  only  slight  effects  of  incoming 
velocities  on  droplet  lifetime  ratios.  The  maximum  difference,  which  occurs  for  the 
leading  droplet  between  the  droplet  lifetime  ratios  at  two  velocities  at  a  H/R  ratio  of  4,  is 
less  than  0.1.  The  same  situation  applies  to  the  cases  at  a  pressure  of  400  atm,  as  shown 
in  Fig.  6.15. 

At  a  higher  pressure  of  400atm,  droplet  interaction  effect  becomes  much  weaker, 
as  already  discussed  above.  In  Fig.  6.15,  the  droplet  lifetime  ratios  of  the  leading  droplets 
are  all  above  0.9  because  of  negligible  interaction  effect.  For  the  trailing  droplet,  the 
droplet  lifetime  ratios  at  a  H/R  ratio  of  8  becomes  larger  than  0.8,  which  is  significantly 
higher  than  its  counterparts  at  a  pressure  of  100  atm  (e.g.,  ~  0.5). 
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Fig.  6.15.  Droplet  lifetime  ratios  vs.  H/R  ratios  at  p=400  atm 
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6.3.5  Droplet  Trajectories 

For  the  droplet  vaporization  case,  the  total  drag  force  acting  on  a  droplet  consists 
of  three  components;  the  friction  drag,  the  pressure  drag,  and  the  thrust  drag.  The  thrust 
drag  results  from  the  surface  blowing.  With  two  droplets  moving  in  tandem,  the  moving 
trajectories  for  the  leading  and  trailing  droplets  are  significantly  different,  especially  at 
low  H/R  ratios.  Fig  6.16  presents  the  droplet  moving  trajectories  at  the  same  pressure  and 
incoming  velocity  of  100  am  and  20  m/s  but  different  H/R  ratios.  The  presence  of  the 
trailing  droplet  performs  only  a  minor  effect  on  the  motion  of  the  leading  droplet,  whose 
trajectory  generally  follows  that  of  an  isolated  droplet.  The  situation  for  the  trailing 
droplet  at  a  H/R  ratio  of  4  is  dramatically  different.  The  trajectory  of  the  trailing  droplet 
deviates  from  that  of  the  leading  droplet  after  a  short  time.  After  that,  the  droplet  moves 
extremely  slowly.  From  50  ns  on,  the  droplet  even  begins  a  slightly  forward  movement, 
corresponding  to  a  forward  drag  force  dominated  by  the  pressure  drag.  After  reaching 
about  85  when  the  leading  droplet  undergoes  break-up  and  is  blown  over  the  top  of  it, 

the  trailing  droplet  begins  to  move  like  an  isolated  droplet  as  expected. 

At  a  H/R  ratio  of  8,  there  is  also  a  transition  time  when  the  movement  of  the 
trailing  droplet  tremendously  slows  down,  but  there  is  no  forward  movement  observed. 
After  the  transition  time,  the  leading  droplet  disappears,  resulting  in  a  lone  second  droplet 
moving  at  much  faster  speeds. 


Time  (|lis) 

Fig.  6.16.  Droplet  trajectories  for  both  the  leading  and  trailing  droplets  at  p=100  atm,  and 
u=20  m/s. 


Time  (|is) 

Fig.  6.17.  Droplet  trajectories  for  both  the  leading  and  trailing  droplets  at  p=400  atm,  and 
u=5  m/s. 
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Pressure  effect  on  droplet  trajectories  are  illustrated  in  Fig.  6.17,  which  presents 
droplet  movements  at  a  pressure  and  an  incoming  velocity  of  400  atm  and  5  m/s 
( Re  ~  120 ),  respectively.  In  Fig.  6.17,  the  only  difference  between  the  droplet  trajectories 
of  the  leading  and  trailing  droplets  occurs  at  a  H/R  ratio  of  4.  This  result  further  proves 
the  conclusion  that  pressure  renders  weak  droplet  interactions.  At  a  H/R  ratio  of  4,  the 
transition  time  for  the  trailing  droplet  begins  at  about  100  jds ,  after  which  point  the 
movement  of  the  trailing  droplet  decelerates  through  the  end  of  the  vaporization  process. 
There  is  no  forward  droplet  movement  observed  for  all  the  cases  shown  in  Fig.  6.17. 

6.4  Summary  of  Results 

Supercritical  vaporization  of  two  oxygen  droplets  moving  in  tandem  in  convective 
hydrogen  environments  is  numerically  investigated.  The  influences  of  the  immediate 
neighboring  droplet  on  droplet  lifetimes,  deformations,  and  trajectories  are  studied  in 
detail. 

Based  on  the  definition  of  the  droplet  surface  at  an  isothermal  line  of  oxygen 
critical  temperature  (154.6K),  there  is  no  droplet  collision  observed  for  all  the  cases 
studied  herein.  However,  a  forward  bag  break-up  of  the  leading  droplet  is  found  when  the 
two  droplets  are  initially  positioned  closely  at  a  H/R  ratio  of  4  and  a  pressure  of  100  atm. 
According  to  mass  fraction  fields,  droplet  collisions  do  happen  at  two  initial  H/R  ratios  of 
4  and  8. 

Deformations  of  the  two  droplets  display  distinct  characteristics.  The  leading 
droplet  confronting  strong  incoming  flow  expands  perpendicular  to  the  flow  direction, 
while  the  trailing  droplet  mainly  deform  along  the  flow  direction. 
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Initially  there  are  two  recirculating  eddies  developed  independently  in  the  wakes 
of  the  two  droplets.  With  an  incoming  velocity  of  10  m/s  and  an  initial  H/R  ratio  of  4,  the 
second  eddy  disappears  before  the  end  of  the  vaporization  process,  but  a  new  eddy  is 
created  behind  the  first  one.  Increasing  the  initial  H/R  ratio  to  8,  however,  changes  the 
flow  fields  dramatically.  In  the  later  case,  the  first  eddy  disappears,  after  the  leading 
droplet  is  completely  evaporated.  The  second  eddy  remains  during  the  entire  vaporization 
process. 

The  existence  of  the  trailing  droplet  produces  only  minor  effects  on  the  lifetimes 
and  trajectories  of  the  leading  droplets.  At  a  pressure  of  100  atm  and  the  initial  H/R  ratios 
of  4  and  8,  the  lifetimes  of  the  trailing  droplet  are  strongly  elongated.  There  is  forward 
movement  observed  for  the  trailing  droplet  at  a  H/R  ratio  of  4  and  a  pressure  of  100  atm. 

Increase  of  pressure  weakens  droplet  interactions,  and  this  conclusion  is 
consistently  verified  in  the  studies  of  temperature  fields,  flow  fields,  droplet  lifetimes, 
and  droplet  trajectories. 
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Fig.  6.3.  Transient  variations  of  temperature  fields  and  droplet  deformatii 
atm,  u=10  m/s,  and  H/R  ratio  of  4  (t=6,  60,  155,  and  200  ps,  respectively). 
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Fig.  6.4.  Transient  variations  of  temperature  fields  and  droplet  deformations  at  p=100 
atm,  u=10  m/s,  and  H/R  ratio  of  8  (t=15,  145,  190,  and  225  ps,  respectively). 
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Fig.  6.5.  Transient  variations  of  temperature  fields  and  droplet  deformations  at  p=100 
atm,  u=10  m/s,  and  H/R  ratio  of  12  (t=20,  60, 100,  and  150  ps,  respectively). 
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Fig.  6.6.  Transient  variations  of  oxygen  mass  fraction  fields  at  p=100  atm,  u=10  m/s,  and 
H/R  ratio  of  4  (t=6,  60, 125,  and  200  ps,  respectively). 
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Fig.  6.8.  Transient  variations  of  flow  fields  at  p=100  atm,  u=20m/s,  and  H/R  ratio  of  4 
(t=5, 40,  65,  85  ps,  respectively). 
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Fig.  6.9.  Transient  variations  of  flow  fields  at  p=100  atm,  u=20m/s,  and  H/R  ratio  of  8 
(t=8, 40,  90,  1 10  ps,  respectively). 
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Fig.  6.10.  Transient  variations  of  flow  fields  at  p=100  atm,  u=20m/s,  and  H/R  ratio  of  12 
(t=15, 40,  60,  80  ps,  respectively). 
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Chapter  7 


CONCLUSIONS  AND  RECOMMENDATIONS 

A  systematic  numerical  investigation  of  supercritical  droplet  vaporization  and 
droplet  interactions  has  been  conducted  based  on  the  transient  conservation  equations  of 
mass,  momentum,  energy,  and  species.  A  unified  treatment  of  real  fluid  behavior  was 
developed,  which  is  capable  of  treating  fluid  states  over  the  entire  thermodynamic 
regime,  ranging  from  dense  liquid  to  dilute  gas.  Fundamental  thermodynamic  theories 
were  applied  for  the  derivations  of  thermodynamic  relationships,  which  naturally  take 
pressure  effects  and  thermodynamic  non-ideality  into  account.  A  preconditioning  scheme 
with  the  dual  time-stepping  integration  technique  was  further  used  to  make  the  numerical 
algorithm  capable  of  solving  small  Mach-number  fluid  flows.  In  the  current  research 
work,  a  modified  Soave-Redlich-Kwong  (SRK)  equation  of  state  was  utilized  for 
evaluating  all  the  thermodynamic  correlations,  which  was  then  incorporated  into  the 
numerical  scheme  to  enhance  its  efficiency  and  robustness.  The  state-of-the-art 
techniques  were  employed  to  estimate  the  transport  properties,  and  full  account  was  taken 
for  the  effect  of  transport  anomaly  at  the  transcritical  region.  The  extended  corresponding 
states  one-fluid  model  combined  with  the  Benedict-Webb-Rubin  (BWR)  equation  of  state 
was  implemented  to  predict  the  viscosity  and  thermal  conductivity  of  the  mixtures,  while 
the  pressure  effect  on  binary  mass  diffusivity  was  corrected  by  the  Takahashi  method. 

A  linear  correlation  relating  mass  vaporization  rate  with  the  difference  of 
chemical  potentials  in  both  liquid  and  gaseous  phases,  which  is  based  on  Onsager’s 
imeversible  thermodynamic  theory,  was  incorporated  into  the  unified  preconditioning 


numerical  algorithm  to  study  droplet  vaporization  in  a  quiescent  environment  in  Chapter 
3.  The  oxygen-hydrogen  system  was  examined  under  both  sub-  and  super-critical 
conditions.  Droplet  cluster  behavior  in  the  inner  region  of  a  droplet  cluster  was  explored 
by  assuming  an  isolated  and  impermeable  bubble  surrounding  each  droplet.  Studies  were 
focused  on  the  pressure  and  temperature  effects  on  droplet  interactions.  Results 
demonstrate  that  pressure  has  a  strong  effect  on  droplet  interactions,  but  the  temperature 
effect  is  minor,  especially  at  high  pressures.  With  the  increase  of  pressure,  the  effect  of 
droplet  interactions  on  droplet  lifetime  decreases.  The  effects  of  pressures  and 
temperatures  on  hydrogen  density  play  a  decisive  role  in  determining  droplet  interactions 
by  ultimately  influencing  the  temperature  and  mass  fraction  gradients  at  the  droplet 
surface.  At  the  later  stage  of  the  vaporization  process,  the  lower  gradient  of  oxygen  mass 
fraction  at  the  droplet  surface  causes  more  oxygen  accumulation,  which  results  in  lower 
thermal  conductivity  of  the  mixture  and  consequently  slower  heat  transfer  rate  into  the 
liquid  droplet.  This  factor  dominates  at  the  later  stage  of  the  droplet  vaporization  process. 

An  oxygen  droplet  vaporization  in  forced-convection  hydrogen  environments  was 
studied  in  Chapter  5,  where  the  droplet  was  assumed  to  reach  the  critical  mixing  state 
immediately  after  being  injected  into  the  combustion  chamber.  The  effects  of  the 
Reynolds  number  and  pressure  on  droplet  deformations,  flow  filed  developments,  droplet 
lifetimes,  and  droplet  dynamics  were  thoroughly  examined.  A  dimensionless  parameter, 

joh ,  which  represents  the  ratio  of  aerodynamic  and  viscous  forces,  was  found  to 
be  the  major  factor  determining  the  droplet  deformations  under  supercritical  conditions. 
Both  Reynolds  number  and  pressure  perform  their  effects  by  changing  this  parameter. 
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At  a  given  pressure,  a  droplet  evaporates  more  rapidly  at  a  higher  Reynolds 
number,  while  the  droplet  lifetime  is  significantly  shortened  at  a  lower  pressure  and  a 
given  Reynolds  number.  A  relationship  was  established  to  accurately  calculate  droplet 
lifetimes. 
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A  linear  relationship  was  reasonably  assumed  for  droplet  velocity,  and  the  linear 
velocity  coefficient  b  was  found  to  proportional  to  Re^'^^  .  Droplet  positions 

calculated  from  this  linear  relationship  agree  well  with  those  directly  from  numerical 
computations. 

b  =  (15<Re<240,2<F/j  <8)  (7.3) 

The  drag  coefficient  in  the  supercritical  vaporization  case  always  decreases  with 
f 

fYl 

time  ( ~  —  )  as  a  result  of  the  weak  Reynolds  number  reduction.  At  a  given 

roj 


Reynolds  number,  the  drag  coefficient  decreases  slightly  with  the  increase  of  pressure. 

Droplet  interactions  in  convective  environments  were  examined  in  Chapter  6, 
where  studies  were  focused  on  two  droplets  moving  in  tandem.  The  goal  is  to  identify  the 
influences  from  the  immediate  neighboring  droplet  on  droplet  deformations,  lifetimes, 
and  trajectories.  Deformations  of  the  two  droplets  display  distinct  characteristics.  The 
leading  droplet  confronting  strong  incoming  flow  expands  perpendicular  to  the  flow 
direction,  while  the  trailing  droplet  deforms  mainly  along  the  flow  direction.  Based  on 
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the  definition  of  the  droplet  surface  at  an  isothermal  line  of  oxygen  critical  temperature 
(154.6K),  there  is  no  droplet  collision  observed  for  all  the  cases  studied  herein.  However, 
a  forward  bag  break-up  of  the  leading  droplet  was  found  when  the  two  droplets  were 
initially  positioned  closely  at  a  ff/R  ratio  of  4  and  a  pressure  of  100  atm.  According  to 
mass  fraction  fields,  droplet  collisions  do  happen  at  two  initial  H/R  ratios  of  4  and  8.  The 
existence  of  the  trailing  droplet  produces  only  slight  effects  on  the  lifetimes  and 
trajectories  of  the  leading  droplets,  while  the  lifetimes  of  the  trailing  droplet  are  strongly 
elongated.  Forward  movement  was  observed  for  the  trailing  droplet  at  a  H/R  ratio  of  4 
and  a  pressure  of  lOOatm.  Increase  in  pressure  weakens  droplet  interactions,  and  this 
conclusion  was  consistently  verified  in  the  studies  of  temperature  fields,  flow  fields, 
droplet  lifetimes,  and  droplet  trajectories. 

Results  in  this  research  work  enhance  the  fundamental  physical  understanding  of 
droplet  vaporization  and  cluster  behavior  at  both  sub-  and  super-critical  conditions.  The 
correlations  regarding  droplet  lifetimes  and  droplet  dynamics  can  be  applied  directly  into 
the  spray  combustion  modeling. 

In  the  current  research  work,  combustion  phenomenon  was  not  included  due  to 
computational  time.  In  the  future  studies,  the  numerical  code  is  expected  to  be 
parallelized  to  take  advantage  of  the  most  advanced  computational  technique,  which 
renders  detailed  studies  of  combustion  and  flame  around  droplets  feasible. 

With  the  increased  computational  efficiency,  this  computational  algorithm  can  be 
further  extended  for  three-dimensional  calculations.  Droplet  interactions  with  two 
droplets  moving  side  by  side  can  be  studied  in  detail  to  complement  our  basic 
understanding  of  this  most  common  droplet  interacting  phenomenon.  In  addition,  the 
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computational  algorithm  can  also  be  used  to  study  three-dimensional  supercritical  spray 
and  mixing  layer,  which  display  features  distinct  from  the  low-pressure  counteiparts 
based  on  the  current  limited  experimental  results. 

The  greatest  challenge  is  to  study  hydrocarbon  fuel  droplet  vaporization  in 
supercritical  convective  environments,  since  the  droplet  is  expected  to  spend  sufficiently 
long  time  in  subcritical  regions.  Here  a  high-pressure  phase-changing  problem  with 
strong  surface  deformation  and  possible  droplet  break-up  has  to  be  treated.  The  internal 
phase-changing  boundary  conditions  are  too  complicated  for  the  current  surface  tracking 
numerical  schemes  to  handle,  as  presented  in  Chapter  2.  The  success  in  this  area  will 
depend  mainly  on  the  future  researches  on  phase-changing  physics  and  the  development 
of  the  simple  and  reliable  phase-changing  correlations. 
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Appendix  A 


THERMODYNAMIC  RELATIONSHIPS 


Based  on  the  Soave-Redlich-Kwong  equation  of  state,  the  following  derivative 
expressions,  which  are  used  extensively  in  the  following  thermodynamic  relationships, 
can  be  directly  derived  as 
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where  the  indexes  i,  j  =  I--- N ,  and  the  derivative  —  {a a)  is  given  in  Appendix  B. 

dT 

It  should  be  noted  that,  according  to  thermodynamics,  these  derivatives  could  be 
derived  based  on  any  equation  of  state. 

The  partial  density  internal  energy  of  species  i  will  then  be  derived.  We  first 

need  to  find  the  expression  for  the  internal  energy  e .  From  the  fundamental 
thermodynamic  theory,  we  have 
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where  the  subscript  0  indicates  a  reference  ideal  state  at  a  low  pressure. 

Utilizing  SRK  equation  of  state  and  the  partial  derivative  relation  (Al),  Eq.  (A4) 
is  integrated,  which  leads  to  the  following  relationship: 
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where  the  partial  derivative  -  is  also  presented  in  Appendix  B. 

oT 

According  to  the  definition  for  the  partial  density  property,  the  partial  density 
internal  energy  ej  can  be  expressed  as 
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In  Addition,  utilizing  Eq.  (A6),  the  internal  energy  of  a  mixture  can  be  related  to 
the  partial  density  internal  energy  as 
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Based  on  thermodynamics, 
ph  =  pe  +  p 


(A8) 


Following  the  definition  for  the  partial  density  property,  the  following  expression 
can  be  found  by  taking  derivative  of  the  partial  density  of  species  i  to  both  sides  of  the 
Eq.  (A8),  and  keeping  temperature  and  all  the  other  partial  densities  constant: 
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It  is  recognized  that  Eq.  (A9a)  is  equivalent  to 
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The  following  thermodynamic  expression  exists,  which  relates  a  partial  density 
property  to  a  partial  mass  property: 
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where  the  parameter  <p  refers  to  any  proper  intensive  thermodynamic  property,  such  as 
enthalpy,  and  internal  energy. 

Substituting  Eq.  (A9c)  into  Eq.  (A9b),  and  taking  use  of  the  fundamental  enthalpy 
expression,  which  can  be  found  in  any  thermodynamics  textbook,  the  following  relation 
concerning  the  partial  mass  enthalpy  /?/  can  be  further  derived: 
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Next,  we  begin  to  derive  the  expressions  for  the  constant  volume  and  the  constant 
pressure  heat  capacities  based  on  the  SRK  equation  of  state. 

The  definition  of  constant  volume  heat  capacity  is 
^de'^ 
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Utilizing  Eq.  (A5),  it  is  straightforward  to  find 
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where  the  derivative  — is  given  in  Appendix  B. 
dT^ 

Following  fundamental  thermodynamic  relationships,  the  constant-pressure  heat 
capacity  can  be  expressed  as 
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In  order  to  find  the  thermodynamic  relationships  regarding  chemical  potential,  the 
partial  density  and  partial  mass  entropy  have  to  be  derived  first. 
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Based  on  the  definition  of  partial  density  entropy,  it  is  found 
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The  partial  mass  entropy  can  be  further  related  to  the  partial  density  entropy  as 
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where  the  partial  mass  volume  is 
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The  chemical  potential  of  species  i  can  be  calculated  as 
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Finally,  the  partial  derivatives  regarding  chemical  potential  can  be  expressed  as 
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Appendix  B 

DERIVATIVE  EXPRESSIONS  IN  SOAVE-REDLICH-KWONG 

EQUATION  OF  STATE 

In  the  Soave-Redlich-Kwong  (SRK)  equation  of  state,  the  terms  aa  and  qm.j  are 
functions  of  temperature.  The  derivatives  of  these  terms  to  temperature  are  derived  as 
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The  second  derivatives  of  parameter  aa  to  temperature  is 
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The  variable  o;,-  for  species  H2,<7://2  >  treated  differently  since  hydrogen  is  a 
quantum  gas.  The  derivatives  of  this  variable  are 
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NUMERICAL  JACOBIANS  IN  PRECONDITIONING  SCHEME 


Jacobian  T  is  derived  as 
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Jacobian  B  is 
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Jacobian  is  expressed  as 
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In  Jacobian  D,  only  those  terms  related  to  axisymmetric  case  are  given.  The 
combustion  terms  can  be  found  in  Shuen  et  al.  (1990). 
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Appendix  D 

NUMERICAL  MODEL  VALIDATION 

Dl.  Validation  of  Equation  of  State 

In  this  research  work,  the  modified  Soave-Redlich-Kwong  (SRK)  equation  of 
state  is  used  for  calculating  all  the  thermodynamic  relationships,  and  the  Benedict-Webb- 
Rubin  (BWR)  equation  of  state  is  used  for  evaluating  the  transport  properties,  including 
thermal  conductivity  and  viscosity.  Figure  Dl  presents  the  relative  errors  associated  with 
these  equations  of  state  for  calculating  the  density  of  oxygen.  The  results  from  the  BWR 
equation  of  state  show  excellent  agreement  with  experimental  data.  Since  the  coefficients 
of  BWR  equation  of  state  are  currently  only  available  for  several  pure  substances,  the 
SRK  equation  of  state  is  used  for  calculating  thermodynamic  relationships.  Figure  Dl 
indicates  that  the  relative  error  associated  with  SRK  equation  of  state  is  only  about  10%, 
which  is  better  than  that  of  the  Peng-Robinson  (PR)  equation  of  state. 
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Figure  Dl  Relative  errors  associated  with  SRK  and  BWR  equations  of  state  for 
calculating  the  density  of  oxygen.  From  Yang  (2000).  Used  with  permission. 
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D2.  Calculation  of  Thermodynamic  Phase  Equilibrium 

Figure  D2  presents  the  thermodynamic  states  of  the  n-pentane  and  carbon  dioxide 
system  at  phase  equilibrium,  which  are  calculated  based  on  the  modified  Soave-Redlich- 
Kwong  equation  of  state.  The  results  are  compared  with  the  available  experimental  data 
(Poettmann  and  Katz  1945)  at  two  temperatures  of  366.7  K  and  394.4  K.  Figure  D2 
clearly  illustrates  the  pressure-dependent  dissolution  of  gaseous  carbon  dioxide  into 
liquid  n-pentane.  The  computational  results  present  excellent  agreement  with 
experimental  data  for  the  gas-phase  equilibrium  compositions,  but  deviate  slightly  from 
the  liquid-phase  experimental  data.  This  is  mainly  caused  by  the  lack  of  appropriate 
experimental  data  for  the  binary  interaction  coefficients  in  the  SRK  equation  of  state.  It 
may  also  be  attributed  to  the  relative  errors  of  the  SRK  equation  of  state  for  estimating 
liquid  density. 


Figure  D2  Comparison  of  computational  and  experimental  results  for  high-pressure  phase 
equilibrium,  n-pentane/carbon  dioxide  system. 
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D3.  Validation  of  Numerical  Model 

The  numerical  model  of  single  droplet  vaporization  in  a  quiescent  environment, 
which  is  developed  in  Chapter  2,  is  validated  here  by  comparing  numerical  results  with 
experimental  data.  Since  an  assumption  of  non-gravity  is  used  in  the  development  of  the 
numerical  model  to  avoid  natural-convection  effect,  only  the  experimental  data  obtained 
under  microgravity  conditions  can  be  used  for  comparison.  The  most  recent  experimental 
results  of  Nomura  et  al.  (1996)  are  chosen  for  model  validation.  During  the  experiment, 
one  n-heptane  droplet  was  suspended  in  a  nitrogen  environment.  The  initial  droplet 
diameter  was  0.6-0.8  mm,  and  pressure  ranged  from  0.1  to  5.0  Mpa.  The  microgravity 
conditions  were  achieved  using  5-m  and  110-m  drop  tower  and  parabolic  flights,  which 
reduced  the  gravity  level  to  the  order  of  10'"-10‘^g.  The  temporal  variation  of  droplet 
diameter  was  measured  with  a  computer-aided  image  analysis  system.  Before  we 
compare  the  numerical  and  experimental  results,  the  following  factors,  which  may  cause 
discrepancies  between  the  results,  should  be  emphasized.  First,  there  were  inevitable 
experimental  errors  associated  with  the  image  analysis  system,  which  were  not  reported. 
Second,  the  experiments  were  conducted  for  a  suspended  droplet,  but  the  numerical 
computation  is  conducted  for  a  floating  droplet.  Recent  estimates  of  Morin  (1999) 
indicate  that  the  effect  of  suspension  fiber  on  the  droplet  vaporization  rate  increases  with 
the  increase  of  pressure  and  temperature,  and  with  the  decrease  of  droplet  size.  Finally, 
the  results  of  Vielli  et  al.  (1996)  demonstrate  that  the  microgravity  conditions  obtained  in 
the  drop  tower  and  parabolic  flight  experiments  can  not  completely  avoid  the  buoyancy 
effect  on  droplet  vaporization  rate,  which  increases  with  the  increase  of  pressure. 
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Figure  D3  presents  the  temporal  variations  of  the  dimensionless  droplet  diameter 
at  a  pressure  of  0.5  Mpa,  and  two  ambient  temperatures  of  660  and  750  K,  respectively. 
Reasonably  good  agreement  has  been  reached  for  both  cases.  At  the  beginning  of  the 
vaporization  process,  both  the  experimental  and  computational  results  illustrate  droplet 
expansion,  which  is  caused  mainly  by  sudden  heating  of  the  droplet.  Figure  D4  compares 
the  numerical  and  experimental  results  at  a  pressure  of  2  Mpa  and  an  ambient 
temperature  of  660  K.  Excellent  agreement  is  reached  in  the  early  stage  of  the 
vaporization  process,  but  the  difference  grows  during  the  vaporization  process.  This 
phenomenon  is  consistent  with  the  conclusions  of  Morin  (1999)  and  Vielli  et  al.  (1996). 
The  reasonably  good  agreement  of  the  droplet  vaporization  rate  between  the  numerical 
and  experimental  results  validates  the  numerical  model  developed  in  this  research  work. 
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Fig.  D3  Comparisons  of  numerical  and  experimental  results  of  n-heptane  droplet 
vaporization  in  nitrogen  environments  at  pressure  of  0.5  Mpa  and  two  ambient 
temperatures  of  660  and  750  K. 


Fig.  D4  Comparison  of  numerical  and  experimental  results  of  n-heptane  droplet 
vaporization  in  nitrogen  environment  at  pressure  of  2  Mpa  and  ambient  temperature  of 
660  K. 


